{'fraction_non_alphanumeric': 0.03966518044923025, 'fraction_numerical': 0.014259274838058384, 'mean_word_length': 4.41992476917816, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimizing squeezing in a coherent quantum feedback network of optical parametric oscillators

Advances in the emerging field of coherent quantum feedback control (CQFC) have led to the development of new capabilities in the areas of quantum control and quantum engineering, with a particular impact on the theory and applications of quantum optical networks. We consider a CQFC network consisting of two coupled optical parametric oscillators (OPOs) and study the squeezing spectrum of its output field. The performance of this network as a squeezed-light source with desired spectral characteristics is optimized by searching over the space of model parameters with experimentally motivated bounds. We use the QNET package to model the network's dynamics and the PyGMO package of global optimization algorithms to maximize the degree of squeezing at a selected sideband frequency or the average degree of squeezing over a selected bandwidth. The use of global search methods is critical for identifying the best possible performance of the CQFC network, especially for squeezing at higher-frequency sidebands and higher bandwidths. The results demonstrate that the CQFC network of two coupled OPOs makes it possible to vary the squeezing spectrum, effectively utilize the available pump power, and overall significantly outperform a single OPO. Additionally, the Hessian eigenvalue analysis shows that the squeezing generation performance of the optimally operated CQFC network is robust to small variations of phase parameters.


I. INTRODUCTION
Feedback control is ubiquitous in classical engineering. However, its extension to the quantum realm has been challenging due to the unique character of the quantum measurement, which requires coupling of the observed quantum system to a classical measurement apparatus. Consequently, measurement-based quantum control has to deal with the fundamental effect of stochastic measurement back action on the quantum system, along with the need to amplify quantum signals up to macroscopic levels and high latency of classical controllers in comparison to typical quantum dynamic time scales [1,2]. An alternative approach that has attracted significant interest in the last decade is coherent quantum feedback control (CQFC) [3][4][5], which considers networks where the quantum system of interest (called the plant) is controlled via coupling (either direct or, more often, through intermediate quantum fields) to an auxiliary quantum system (called the controller). CQFC schemes utilize coherent quantum signals circulating between the plant and controller, thus avoiding the need for signal amplification and associated excess noise. Also, both plant and controller can evolve on the same time scale, which eliminates the latency issues. Due to these advantages, CQFC makes it possible to engineer quantum networks with new and unique characteristics [4][5][6][7].
The theoretical foundation of CQFC is a powerful framework based on input-output theory, which is used for modeling networks of open quantum systems connected by electromagnetic fields [8][9][10][11] (see also [3][4][5] for reviews). Moreover, recent developments, including the SLH formalism [12][13][14], the quantum hardware description language (QHDL) [15], and the QNET software package [16], have added important capabilities for, respectively, modular analysis, specification, and simulation of such quantum optical networks. Together, the existing theoretical tools enable efficient and automated design and modeling of CQFC networks.
Proposed and experimentally demonstrated applications of CQFC include the development of autonomous devices for preparation, manipulation, and stabilization of quantum states [17][18][19][20][21][22], disturbance rejection by a dynamic compensator [23], linear-optics implementation of a modular quantum memory [24], generation of optical squeezing [25][26][27][28], generation of quantum entanglement between optical field modes [29][30][31][32][33][34], coherent estimation of open quantum systems [35,36], and ultra-low-power optical processing elements for optical switching [37][38][39] and analog computing [40,41]. In addition to tabletop bulk-optics implementations, CQFC networks have been also implemented using integrated silicon photonics [42] and superconducting microwave devices [43,44]. Squeezed states of light [45][46][47][48] have found numerous applications in quantum metrology and quantum information sciences, including interferometric detection of gravitational waves [49,50], continuous-variable quantum key distribution (CV-QKD) [51][52][53][54][55][56], generation of Gaussian entanglement [55][56][57], and quantum computing with continuousvariable cluster states [58][59][60][61]. Different applications require squeezed states with different properties. For example, detectable gravitational waves are expected to have frequencies in the range from 10 Hz to 10 kHz, and, consequently, quadrature squeezed states used to increase the measurement sensitivity in interferometric detectors should have a high degree of squeezing at sideband frequencies in this range. On the other hand, in CV-QKD the secure key rate is proportional to the bandwidth of squeezing, and hence it would be useful to generate states with squeezing bandwidth extending to 100 MHz or even higher. It would be also of interest to extend the maximum of squeezing to high sideband frequencies.
In recent years, there have been remarkable advances in the generation of squeezed states [62][63][64][65][66][67][68][69][70][71][72][73][74], however, achieving significant control over the squeezing spectrum still remains an ongoing effort. In 2009, Gough and Wildfeuer [25] proposed to enhance squeezing in the output field of a degenerate optical parametric oscillator (OPO) by incorporating the OPO into a CQFC network, where a part of the output beam is split off and then fed back into the OPO. Iida et al. [26] reported an experimental demonstration of this scheme, while Német and Parkins [28] proposed to modify it by including a time delay into the feedback loop. Another significant modification of this scheme was proposed and experimentally demonstrated by Crisafulli et al. [27], who included a second OPO to act as the controller, with the plant OPO and the controller OPO coupled by two fields propagating between them in opposite directions. Due to the presence of quantum-limited gains in both arms of the feedback loop, this CQFC network has a very rich dynamics. In particular, by tuning the network's parameters it is possible to significantly vary the squeezing spectrum of its output field, for example, shift the maximum of squeezing from the resonance to a high-frequency sideband [27].
The full range of performance of the CQFC network of two coupled OPOs as a squeezed-light source, however, still remains to be explored. In this paper, we study the limits of the network's performance by performing two types of optimizations: (1) maximizing the degree of squeezing at a chosen sideband frequency and (2) maximizing the average degree of squeezing over a chosen bandwidth; in both cases, the searches are executed over the space of network parameters with experimentally motivated bounds. To maximize the chances of finding a globally optimal solution, we use the PyGMO package of global optimization algorithms [75] and employ a hybrid strategy which executes in parallel eight searches (using seven different global algorithms). Before each optimization is completed, the searches are repeated multiple times, and intermediate solutions are exchanged between them after each repetition. This strategy enabled us to discover that the CQFC network, when optimally operated, is capable of achieving a remarkably high degree of squeezing at sideband frequencies and bandwidths as high as 100 MHz, with a very effective utilization of the available pump power. We also find that the obtained optimal solutions are quite robust to small variations of phase parameters.

II. BACKGROUND
The derivations in this section largely follow those in Refs. [25,27], with some additional details and modifications.
A. Input-output model of a quantum optical network Consider a network of coupled linear and bilinear optical elements such as mirrors, beam-splitters, phase-shifters, lasers, and degenerate OPOs. The quantum theory of such a network considers quantized cavity field modes which are coupled through cavity mirrors to external (input and output) quantum fields [9][10][11]. Let n be the number of the network's input ports (equal to the number of output ports) and m be the number of cavities (in this model, we assume that each cavity supports one internal field mode). Let a, a in , and a out denote vectors of boson annihilation operators for, respectively, the cavity modes, the input fields, and the output fields: Assuming that all input fields are in the vacuum state, the network is fully described by the (S, L, H) model (also called the SLH model) [12][13][14], which includes the n × n matrix S that describes the scattering of external fields, the n-dimensional vector L that describes the coupling of cavity modes and external fields, and the Hamiltonian H that describes the intracavity dynamics. For the model considered here, elements {S ij } of S are c-numbers, while H and elements {L i } of L are operators on the combined Hilbert space of all cavity modes in the network. The Heisenberg equations of motion (also known as quantum Langevin equations) for the cavity mode operators {a (t)} are ( = 1) Here, L L is the Lindblad superoperator: and Γ l is the noise operator: where a † in = [a † in,1 , . . . , a † in,n ] and L † = [L † 1 , . . . , L † n ] are row vectors of respective Hermitian conjugate operators. The generalized boundary condition for the network is For the type of networks that we consider, elements of L are linear in annihilation operators of the cavity modes, i.e., where K is an n × m complex matrix with elements {K i = [L i , a † ]}, and the Hamiltonian has the bilinear form: where a † = [a † 1 , . . . , a † m ] and a ‡ = a †T are, respectively, row and column vectors of boson creation operators for the cavity modes, Ω is an m × m Hermitian matrix, and W is an m × m complex matrix. With such L and H, the Heisenberg equations of motion (2) take the form: To obtain the transfer-matrix function from input to output fields, we seek the solution of Eq. (8) in the frequency domain. Using the Fourier transform, we define: where b(t) stands for any element of a(t), a in (t), and a out (t).
The field operators are in the interaction frame, and therefore ω is the sideband frequency (relative to the carrier frequency). We also use the double-length column vectors of the form: where b(ω) stands for either of a(ω), a in (ω), and a out (ω).

B. Squeezing spectrum
Consider the quadrature of the ith output field in time and frequency domains: where θ is the homodyne phase. The power spectral density of the quadrature's quantum noise (commonly referred to as the squeezing spectrum) is [45,46]: where : : denotes the normal ordering of boson operators and x, y = xy − x y . Since all input fields are in the vacuum state, X i (ω, θ) = X i (ω , θ) = 0, and one obtains: where By substituting Eqs. (16) into Eqs. (20) and evaluating expectation values for vacuum input fields, one obtains: In this work, we are only concerned with squeezing properties of the field at one of the output ports. We will designate this port as corresponding to i = 1 and denote the squeezing spectrum of this output field as P(ω, θ) = P 1 (ω, θ). In squeezing generation, the figure of merit is the quantum noise change relative to the vacuum level, measured in decibels, and since P vac (ω, θ) = 1, the corresponding spectral quantity is Q(ω, θ) = 10 log 10 P(ω, θ).
Negative values of Q correspond to quantum noise reduction below the vacuum level (i.e., squeezing of the quadrature uncertainty). The maximum degree of squeezing corresponds to the minimum value of Q. The maximum and minimum of P(ω, θ) as a function of θ, are power spectral densities of the quantum noise in antisqueezed and squeezed quadrature, respectively. Analogously to Eq. (22), logarithmic spectral measures of antisqueezing and squeezing for the two quadratures are defined as Q ± (ω) = 10 log 10 P ± (ω), respectively. Expressing M(ω) as M(ω) = |M(ω)|e iθ M (ω) and using Eq. (19), it is easy to find (we omit the subscript i = 1 for simplicity): with anti-squeezed and squeezed quadrature corresponding to θ = θ M (ω)/2 and θ = [θ M (ω) − π]/2, respectively. Note that, in general, these optimum values of the homodyne phase θ depend on the sideband frequency ω, so, for example, if the goal is to maximize the degree of squeezing at a particular sideband frequency ω opt , then the optimum phase value θ opt = [θ M (ω opt ) − π]/2 should be selected accordingly.

III. SQUEEZING FROM A SINGLE OPO
A network that produces squeezed light by means of a single degenerate OPO [47] is schematically shown in Fig. 1. The OPO consists of a nonlinear crystal enclosed in a Fabry-Pérot cavity. The pump field for the OPO is assumed to be classical and not shown in the scheme. Each partially transparent mirror in the network (including cavity mirrors and a beamsplitter) has two input ports and two output ports. A vacuum field enters into each input port. The OPO cavity has a fictitious third mirror to model intracavity losses (mainly due to absorption in the crystal as well as scattering and Fresnel reflection at the crystal's facets). The beamsplitter B models losses in the output transmission line (e.g., due to coupling into a fiber) and inefficiencies in the homodyne detector (not shown) used to measure the squeezing spectrum of the output field. Taking into account all optical elements, the network is modeled as having four input ports, four output ports, and one cavity mode (n = 4, m = 1). Parameters of the single OPO network are described in Table I. With ξ = |ξ|e iθ ξ , there is a total of seven real parameters. Note that we use angular frequencies throughout this paper. For each cavity mirror, the leakage rate is where T i is the power transmittance of the ith mirror, c is the speed of light, and l eff is the effective cavity length (taking into account the length and refractive index of the crystal). To simplify the notation, we also use alternative parameters: to denote the total leakage rate (including losses) from the cavity, and to denote, respectively, the transmittivity and reflectivity of the beamsplitter. The QNET package [16] is used to derive the (S, L, H) model of the network, and the resulting components of the model are where a is the annihilation operator of the cavity field mode. Using the formalism of Sec. II A, we obtain: Ω = ω 0 , W = ξ, where we defined auxiliary parameters: These results make it straightforward to analytically compute the transfer-matrix functionZ(ω) of Eq. (14). Since we are only interested in squeezing properties of the field at the output port 1, it is sufficient to use only the respective rows of matrices Z − (ω) and Z + (ω), i.e., By substituting elements of Z − 1 (ω) and Z + 1 (ω) into Eqs. (21), we obtain: where T B = t 2 B is the power transmittance of the beam splitter. Using Eq. (19), the resulting squeezing spectrum is The spectra for anti-squeezed and squeezed quadrature are obtained as the maximum and minimum (cf. Eq. (23)) of P(ω, θ) in Eq. (30) for ϕ = tan −1 [γω 0 /µ(ω)] and ϕ = tan −1 [γω 0 /µ(ω)] + π, respectively, and are given by In order to compare the theoretical spectra with experimental data, it is common to express the pump amplitude as where P is the OPO pump power and P th is its threshold value. Analogously to the scaled pump amplitude x = 2|ξ|/γ, it is convenient to use scaled frequencies Ω = 2ω/γ and Ω 0 = 2ω 0 /γ. With this notation, Eq. (31) takes the form: where ρ = κ 2 /γ = T 2 /(T 1 + T 2 + L) is the escape efficiency of the cavity, L = T 3 denotes the intracavity power loss, and y 2 = x 2 + Ω 2 − Ω 2 0 . In the case of zero detuning, ω 0 = 0, the squeezing spectrum of Eq. (30) becomes The corresponding spectra for anti-squeezed and squeezed quadrature are obtained for ϕ = 0 and ϕ = π, respectively. They can be expressed by taking Ω 0 = 0 in Eq. (33), which reproduces the familiar result [46,47]: The spectra of Eq. (35) have Lorentzian shapes with maximum (for anti-squeezing) and minimum (for squeezing) at the resonance (zero sideband frequency), and with the degree of squeezing rapidly decreasing as the sideband frequency increases. For applications such as CV-QKD, it would be valuable to significantly extend the squeezing bandwidth. It would Noise power (dB) Anti-squeezing Squeezing Squeezing spectra of the output light field from a single OPO network with different values of the cavity's frequency detuning ω0/2π (given in the legend). Logarithmic power spectral densities of the quantum noise in anti-squeezed and squeezed quadrature, Q ± (ω) = 10 log 10 P ± (ω), are shown versus the sideband frequency ω/2π for P ± (ω) of Eq. (33). The values of network parameters are listed in the text.
be also of interest to achieve a maximum degree of squeezing (i.e., a minimum value of P − ) at a high-frequency sideband. Therefore, we investigate whether such modifications of the squeezing spectrum are possible by using a nonzero value of the cavity's frequency detuning. Consider a single OPO with a set of experimentally motivated parameters: pump power P = 1.5 W, pump wavelength λ p = 775 nm, and signal wavelength λ s = 1550 nm; an MgO:PPLN crystal with length l c = 20 mm, refractive index (at λ s ) n s = 2.1, and effective nonlinear coefficient d eff = 14 pm/V; a Fabry-Pérot cavity with effective length l eff = 87 mm, left mirror reflectance R 1 = 0.98 (T 1 = 0.02, κ 1 /2π ≈ 5.484 MHz), right mirror reflectance R 2 = 0.85 (T 2 = 0.15, κ 2 /2π ≈ 41.132 MHz), intracavity loss L = 0.02 (κ 3 /2π ≈ 5.484 MHz), and total leakage rate γ/2π ≈ 52.1 MHz; output transmission line loss L tl = R B = 0 (T B = 1). These parameters correspond to OPO's threshold power P th ≈ 14.86 W and scaled pump amplitude x = P/P th ≈ 0.318. Using these parameters, we compute the squeezing spectra P ± (ω) of Eq. (33) for three detuning values: ω 0 /2π = {0, 25, 50} MHz. The resulting logarithmic spectra Q ± (ω) = 10 log 10 P ± (ω) for antisqueezed and squeezed quadrature are shown in Fig. 2. These results indicate that, while the use of nonzero detuning can increase the degree of squeezing at higher-frequency sidebands as compared to the case of ω 0 = 0, this increase is very small. Also, no improvement in the squeezing bandwidth (quantified as the average degree of squeezing over a selected bandwidth) is achieved through the use of nonzero detuning. These observations motivate us to explore the use of the CQFC network with two coupled OPOs as a light source with the potential to generate a widely tunable squeezing spectrum.   Real Phase shift of the first phase shifter φ2 Real Phase shift of the second phase shifter θ1 Real Rotation angle of the first beamsplitter θ2 Real Rotation angle of the second beamsplitter θ3 Real Rotation angle of the third beamsplitter

IV. SQUEEZING FROM A NETWORK OF TWO COUPLED OPOS
The CQFC network that includes two coupled degenerate OPOs [27] is schematically shown in Fig. 3. Each OPO consists of a nonlinear crystal enclosed in a Fabry-Pérot cavity.
Pump fields for both OPOs are assumed to be classical and not shown in the scheme. From the control theory perspective, OPO1 is considered to be the plant and OPO2 the (quantum) controller. Each partially transparent mirror in the network (including cavity mirrors and beamsplitters) has two input ports and two output ports. A vacuum field enters into each input port, except for two input ports of cavity mirrors used for the feedback loop between the plant and controller. Each OPO cavity has a fictitious third mirror to model intracavity losses. Beamsplitters B1 and B2 represent the light diverted to lock the cavities as well as losses in optical transmission lines between the OPO cavities. Beamsplitter B3 represents losses in the output transmission line (e.g., due to coupling into a fiber) and inefficiencies in the homodyne detector (not shown) used to measure the squeezing spectrum of the output field. Phase shifters P1 and P2 are inserted into transmission lines between the OPOs to manipulate the interference underlying the CQFC control. Taking into account the feedback loop between the plant and controller, the network is modeled as having seven input ports, seven output ports, and two cavity modes (n = 7, m = 2).
Parameters of the network of two coupled OPOs are listed in Table II. With ξ p = |ξ p |e iθp and ξ c = |ξ c |e iθc , there is a total of 17 real parameters. The relationship between leakage rate and power transmittance of a cavity mirror is given, similarly to Eq. (25), by where T pi (T ci ) is the power transmittance of the ith mirror and l p,eff (l c,eff ) is the effective cavity length for the plant (controller). To simplify the notation, we also use alternative parameters: to denote the total leakage rate (including losses) from, respectively, the plant and controller cavities, to denote, respectively, the transmittivity and reflectivity of each beamsplitter, and to denote the total phase shift for the feedback roundtrip path. Similarly to Eq. (32), we also define the scaled pump amplitudes x p and x c for the plant and controller OPOs, respectively: where P p (P c ) is the OPO pump power and P p,th (P c,th ) is its threshold value for the plant (controller). The QNET package [16] is used to derive the (S, L, H) model of the network, and the resulting components of the model are where a p and a c denote, respectively, the annihilation operators of the plant's and controller's cavity field modes, and we defined auxiliary parameters: By comparing Eq. (43) to the corresponding Hamiltonian without feedback: we observe that two main effects induced by feedback are (1) the appearance of an effective interaction between the plant's and controller's cavity modes, governed by the term i 2 ν 12 a † p a c + H.c., and (2) the modification of the plant detuning by Im ν which is proportional to sin φ.
Using the formalism of Sec. II A, we obtain: where we used additional auxiliary parameters: It is possible to analytically invert the matrixȂ + iωI 4 in order to obtain the transfer-matrix functionZ(ω) and squeezing spectrum P(ω, θ) in analytic form. However, the resulting expressions are too complicated and visually uninformative to be shown here. For practical purposes, it is more efficient to numerically evaluateZ(ω) and P(ω, θ) for any given set of parameter values.

V. SQUEEZING OPTIMIZATION PROCEDURE
A. Objective function In order to quantitatively investigate the tunability of the squeezing spectrum in the CQFC network of two coupled OPOs, we numerically optimize the degree of squeezing at various sideband frequencies. Specifically, we minimize the objective function of the form: where ω opt is the selected sideband frequency. The first term in Eq. (45) is the minimum of the squeezing spectrum at ω opt , while the second term is the uncertainty product times the weight parameter g. This second term is included in order to eliminate solutions with a very large uncertainty of the antisqueezed quadrature. In all optimization results shown below, the weight parameter is g = 0.001. With such a small value of g, the difference between the values of J and P − is always insignificant, and therefore, for the sake of simplicity, we refer to the problem of minimizing J as squeezing optimization.
All solutions encountered during a search are checked to satisfy the Routh-Hurwitz stability criterion [76], i.e., that all eigenvalues of the matrixȂ in Eq. (14) have negative real parts. Any unstable solution is eliminated from the consideration by assigning to it a very large objective value (J = 10 6 ).

B. Optimization variables
For a given ω opt , the objective J is a function of the network parameters -seven real parameters for the single OPO network: and 17 real parameters for the CQFC network of two coupled OPOs: Recall that, for the single OPO network, L = T 3 is the intracavity power loss and L tl = R B is the power loss in the output transmission line. Similarly, for the CQFC network of two coupled OPOs, L p = T p3 and L c = T c3 are the intracavity power losses for the plant and controller OPOs, respectively, and L i = r 2 i (i = 1, 2, 3) are power losses in the transmission lines. In cases where the two intracavity loss values are equal, we denote L in = L p = L c , and where the three transmission line loss values are equal, we denote L out = L 1 = L 2 = L 3 .
Numerical simulations demonstrate that an increase in any of the losses always leads to a deterioration of squeezing, and therefore if a loss parameter can vary in a specified interval [L l , L u ], an optimization will always converge to the lower bound L l . Therefore, it makes sense to to exclude the loss parameters from the optimization variables, i.e., to execute each optimization with all loss parameters having preassigned fixed values (of course, these values can vary from one optimization run to another to explore various experimentally relevant regimes). Consequently, there remain five optimization variables for the single OPO network: and 12 optimization variables for the CQFC network of two coupled OPOs: Each optimization variable z can vary in an interval [z l , z u ] (where z l is the lower bound and z u is the upper bound). The bound intervals are

C. Optimization methodology
Preliminary optimization runs using local algorithms (e.g., Sequential Least Squares Programming) demonstrated that different choices of initial parameter values resulted in different solutions of varying quality. These results mean that the fitness landscape contains multiple local optima. In order to reach a solution of very high quality, we decided to use global search methods. Specifically, we used PyGMO, a suite of global (stochastic) algorithms [75]. Since these global   algorithms are heuristic in nature, they do not guarantee the convergence to a global optimum; in fact, as shown in Table III, while multiple global methods are capable of finding high-quality solutions, the performance varies between different algorithms as well as between optimizations with different values of ω opt for the same algorithm.
To maximize the chances of finding a globally optimal solution, we employed a hybrid strategy, where each optimization executes in parallel eight searches (using seven different global algorithms), with a fully connected topology of solution exchanges between them. These eight searches include two instances of Artificial Bee Colony and one instance of each: Covariance Matrix Adaptation Evolution Strategy, Differential Evolution variant 1220, Differential Evolution with p-best crossover, Improved Harmony Search, Particle Swarm Optimization variant 5, and Compass Search guided by Monotonic Basin Hopping. Each optimization uses the population size of N pop = 30 for each of the global searches, and the evolutions are repeated N ev = 30 times (with solution exchanges between the searches after the completion of each evolution except the last one); the algorithm parameters (the number of no improvements before halting the optimization, N stop , the number of generations, N gen , and the number of iterations, N iter ) used in the searches are the same as those shown in Table III for individual algorithms. As indicated by the results in Table III, this hybrid strategy consistently finds the best solution, as compared to any individual algorithm. Multiple trials with larger values of N pop , N ev , N gen , and N iter did not typically result in an improvement of the solution quality, and thus did not warrant the increased run time.

VI. SQUEEZING OPTIMIZATION RESULTS
First of all, we would like to compare the performance of the CQFC network of two coupled OPOs versus that of the single OPO network, in terms of the maximum degree of squeezing achievable under comparable conditions. Figures 4 and 5 show the optimized degree of squeezing, Q − (ω opt ), at ω opt /2π = 100 MHz, for both networks, versus the upper limits on various network parameters (T u and x u in Fig. 4, and ω u and x u in Fig. 5), with constant loss values: L = L in = 0.01, L tl = L out = 0.1. We observe that the CQFC network of two coupled OPOs generates stronger squeezing than the single OPO network, even as total losses in transmission lines in the former are three times larger than those in the latter (30% versus 10%). In both networks, the maximum degree of squeezing increases with both T u (more light is allowed to leave the cavities) and x u (higher pump power), with these increases being roughly linear for the single OPO network and faster than linear in the CQFC network of two coupled OPOs. These results demonstrate that the feedback makes it possible to more effectively utilize the available pump power. Figure 5 also shows that, for both networks, the maximum degree of squeezing is independent of the upper limit ω u on the cavity detuning frequency; furthermore, we found that in most cases the maximum degree of squeezing is actually achieved with zero detuning. In all results shown below, optimizations used the upper limit value ω u /2π = 100 MHz.
We also investigate the dependence of the maximum degree of squeezing, Q − (ω opt ), on the sideband frequency ω opt at which it is optimized. This dependence is shown in Fig. 6, for both networks, for different values of transmission line losses and pump amplitude bound. We observe that the CQFC network of two coupled OPOs not only generates stronger squeezing than the single OPO network, but that the degradation of squeezing associated with the increase of ω opt is substantially slower in the former than in the latter. The capability of the CQFC network to moderate the degradation of squeezing at higher values of ω opt is associated with a rather abrupt change in the regime of network operation, which is manifested by a rapid change in the slope of the curves in subplots (a) and (c) of Fig. 6.
To explore further the emergence of this new operation regime, we focus on the CQFC network of two coupled OPOs, with Fig. 7 showing the dependence of the maximum degree of squeezing on ω opt for more values of transmission line losses. We observe that the value of ω opt at which the operation regime switches, increases with both x u and L out . The difference between the curve slopes in the low-ω opt and highω opt regimes decreases as L out increases.
To understand the physical differences between operations of the CQFC network in the low-ω opt and high-ω opt regimes, we consider the dependence of the optimal values of power transmittances of cavity mirrors, T p1 , T p2 , T c1 , and T c2 , on ω opt . This dependence is shown in Fig. 8 for optimizations with T u = 0.9, L in = 0.01, and various values of L out and x u . First, we see that the optimal values of T p2 and T c1 are constant over the entire range of ω opt values; specifically, T p2 = 0.9 is at the upper bound, which corresponds to the maximum flow from the plant cavity to the 1st output field (the one whose squeezing properties are measured), and T c1 = 0 is at the lower bound, which corresponds to the minimum flow from the controller cavity to the 5th output field (the one which is not used for either squeezing measurement or feedback). In contrast to this simple behavior of the optimal values of T p2 and T c1 , the optimal values of T p1 and T c2 , which regulate the feedback between the plant and controller OPOs, demonstrate much more intricate dependence on ω opt . The optimal value of T p1 and especially that of T c2 undergo a substantial and rather abrupt change at the critical ω opt value at which the network's operation switches between the low-ω opt and high-ω opt regimes. As ω opt increases through the critical point, T p1 changes from a lower to a higher value, while T c2 decreases from the upper bound T c2 = 0.9 to a much lower value. In other words, the low-ω opt optimal regime is characterized by the maximum flow of light from the controller to the plant and a much lower flow in the opposite direction, while the high-ω opt optimal regime is characterized by roughly similar flows of light in both directions. These patterns characterizing the regimes of optimal network operation, their dependencies on pump and loss parameters, and the rapid switch between the regimes, are quite non-intuitive, and finding them would be rather unlikely without the use of a stochastic global search that explores vast areas of the fitness landscape. The sharp change of the optimal value of T c2 associated with the regime switch makes it easy to identify the sideband frequency ω opt , at which the high-ω opt regime commences (the precision of determining the ω opt /2π values is limited by the sampling interval, which is 2 MHz in our data). The values of ω opt /2π are shown in Table IV for T u = 0.9, L in = 0.01, and various values of L out and x u . We see that ω opt increases monotonously with both L out and x u .
We also note that the optimal values of the scaled pump amplitudes, x p and x c , are almost always at (or very close to) the upper bound x u , i.e., in either regime the optimally operated CQFC network usually uses all the pump power it can get. The maximum use of the pump power is also observed for the optimal operation of the single OPO network. Indeed, as seen in Fig. 9, for both networks, the optimal values of the pump power are virtually independent of ω opt and losses, while they scale quadratically with x u . Due to the rapid growth of the optimal pump power with x u , only values x u ≤ 3 should be considered realistic for the optimal operation with a typical tabletop experimental setup considered in this paper.
Next, we investigate the squeezing spectrum Q − (ω) gen-erated under the optimal operation of either network for various values of ω opt , x u , and transmission line losses. Figure 10 shows Q − (ω) for both networks for various values of ω opt , L tl = L out , and x u . We see that the optimally operated single OPO network generates exactly the same Lorentzian squeezing spectrum for any choice of ω opt . In contrast, the CQFC network of two coupled OPOs is capable of generating diverse squeezing spectra, with the specific spectral shape varying to fit the selected value of ω opt , and overall generates much stronger squeezing over a major portion of the spectrum (especially, at frequencies around ω opt ). Interestingly, the capability of the CQFC network to generate a squeezing spectrum Q − (ω) that has the minimum at ω = ω opt is attained only if the selected value of ω opt is within the highω opt regime of optimal network operation, i.e., ω opt ≥ ω opt (for a given set of bound and loss values). Conversely, as seen for ω opt /2π = 5 MHz in all subplots of Fig. 10 and for ω opt /2π = 25 MHz in subplots (h) and (i) of Fig. 10, the squeezing spectrum has the minimum at ω = 0 if ω opt is within the low-ω opt regime of optimal network operation. Finally, we explore further the dependence of the optimized degree of squeezing, Q − (ω opt ), on the intracavity and transmission line losses for the CQFC network of two coupled OPOs. Figure 11 shows Q − (ω opt ) versus (a) L in and L 3 (with L 1 = L 2 = 0.1), and (b) L in and L out . The situation when L 3 = L 1 = L 2 is practically relevant since L 3 includes, in addition to losses in the output transmission line, inefficiencies in the homodyne detector used to measure the squeezing spectrum of the output field. The results shown in Fig. 11 confirm that any increase in losses is detrimental to squeezing and quantify this relationship.

VII. CORRELATIONS BETWEEN OPTIMAL VALUES OF PHASE VARIABLES
Since the phase parameters play a significant role in tuning the quantum interference that governs the CQFC network's performance, an interesting question is whether their optimal values are correlated. Optimal values of a parameter can be cast as a vector each element of which corresponds to a distinct value of ω opt , and correlations can be computed between pairs of such vectors. Specifically, we computed the Pearson correlation coefficient for all six pairs of four phase variables (θ p , θ c , φ 1 , φ 2 ), and found that substantial correlations only exist between φ 1 and φ 2 . Table V shows the values of the Pearson correlation coefficient r(φ 1 , φ 2 ), computed for T u = 0.9, L in = 0.01, and various values of L out and x u . The correlation in Table V generally decreases as L out and x u increase. This trend can be compared to the one observed in Table IV where the sideband frequency ω opt at which the high-ω opt regime commences increases as L out and x u increase. Since the vectors φ 1 (ω opt ) and φ 2 (ω opt ) contain elements corresponding to both operation regimes, the trend observed in Table V implies that the correlation r(φ 1 , φ 2 ) generally decreases as the number of vector components corresponding to the high-ω opt regime decreases. A plausible explanation of this behavior is that the correlation between the two phase variables is higher in the high-ω opt regime. To test this hypothesis, we computed the Pearson correlation coefficient r(φ 1 , φ 2 ) for the pair of vectors φ 1 = φ 1 (ω opt ≥ ω opt ) and φ 2 = φ 2 (ω opt ≥ ω opt ) that include only elements corresponding to the high-ω opt regime. The values of r(φ 1 , φ 2 ) are shown in Table VI for T u = 0.9, L in = 0.01, and various values of L out and x u . The correlations in Table VI are consistently larger than 0.5, and, furthermore, we find that r(sin φ 1 , sin φ 2 ) = 1.0 and r(cos φ 1 , cos φ 2 ) = −1.0 (up to numerical precision) for all considered values of L out and x u . These findings indicate a significant degree of concerted action in how the CQFC network of two coupled OPOs operates in the high-ω opt regime.

VIII. ROBUSTNESS OF OPTIMAL SOLUTIONS
Any practical implementation of a quantum optical network inevitably involves imprecisions and imperfections, which may affect the desired performance. This issue is of especial importance in a CQFC network, which relies on a precise quantum interference between the pump and controller fields to manipulate the properties of the output field (see, for example, the superposition of a p and a c in the first element of the L vector in Eq. (42)). This interference depends on the values of phase variables, and a key question is how robust is an optimal solution to small variations in these values. To analyze this robustness, we computed the Hessian of the objective function J with respect to the phase variables, for a variety of optimal sets of network parameters.
For the single OPO network, J depends on one phase variable θ ξ , and the Hessian H has one element ∂ 2 J/∂θ The numerical analysis shows that the Hessian is zero (up to numerical precision) for all of these optimal solutions. Therefore, small fluctuations in the value of the pump phase θ ξ should have no effect on the optimized degree of squeezing.
For the CQFC network of two coupled OPOs, J depends on four phase variables (θ p , θ c , φ 1 , φ 2 ), and the Hessian H is a 4 × 4 matrix of second-order derivatives. We The numerical analysis shows that two of the Hessian eigenvalues (h 3 and h 4 ) are zero (up to numerical precision) for all of these optimal solutions. Therefore, robustness to small phase variations is determined by two nonzero Hessian eigenvalues (h 1 and h 2 ). Figures 12-14 show these nonzero Hessian eigenvalues as functions of ω opt and x u (Fig. 12), ω opt and L out (Fig. 13), and x u and L out (Fig. 14). We see that h 1 is typically much larger than h 2 , and hence the magnitude of h 1 is the main factor determining the robustness properties of the optimal solutions. The dependence of h 1 on ω opt , seen in Figs. 12 and 13, demonstrates a significant difference in robustness properties between the low-ω opt and high-ω opt regimes. The low-ω opt regime is intrinsically robust for a broad range of parameter values. In the high-ω opt regime, a reasonable degree of robustness is achieved for x u ≥ 0.2 (i.e., for pump powers above 4 W for the OPO parameters considered here). Larger losses in transmission lines (L out ≥ 0.1) also enhance robustness.
The four components of the Hessian eigenvector e 1 (which corresponds to the largest eigenvalue h 1 ) are shown in Fig. 15 versus ω opt . They also exhibit an abrupt change associated with the switch of the optimal operation regime at ω opt . In the low-ω opt regime, the eigenvector component corresponding to φ 2 has the largest value and the rest of the components have smaller absolute values, but none is negligible. In the highω opt regime, the components corresponding to φ 1 and φ 2 have similar values, while the components corresponding to θ p and θ c are close to zero. These results are consistent with the find- ings that the low-ω opt regime is characterized by the maximum flow of light passing through the phase shifter P2 (from the controller to the plant), while the high-ω opt regime is characterized by roughly similar flows of light passing through the phase shifters P1 and P2 (in both directions). The decrease of the optimized degree of squeezing due to small variations of phase parameters can be quantified using the computed Hessian eigenvalues or, alternatively, via direct Monte Carlo averaging over a random distribution of phase variable values. Figure 16 shows the optimized degree of squeezing, Q − (ω opt ), for the CQFC network of two coupled OPOs (with ω opt /2π = 100 MHz, x u = 0.2, T u = 0.9, and various values of L out ), versus the standard deviation of phase uncertainty, σ phase (for simplicity, we assume a normal distribution with zero mean and the same value of σ phase for uncertainty in each of the four phase variables). We see a good agreement between the Hessian-based and Monte Carlo computations for σ phase ≤ 0.1 (and even for σ phase ≤ 0.2 for L out ≥ 0.1). We also see that the deterioration of squeezing induced by phase variations is quite tolerable for σ phase ≤ 0.1 (especially, for L out ≥ 0.1). Note that our squeezing optimization procedure does not explicitly include a robustness requirement, and hence the observed high level of robustness might be surprising, but it is likely related to the natural tendency of stochastic optimization algorithms to eliminate solutions that are very sensitive to small parameter variations.

IX. SQUEEZING BANDWIDTH OPTIMIZATION
In CV-QKD with squeezed states, the secure key rate is proportional to the bandwidth of squeezing. Therefore, we also explored the capability of the CQFC network of two coupled OPOs to generate output states with high squeezing bandwidth, by optimizing the average degree of squeezing over a frequency interval [0, ω B ], for various values of ω B . Specifically, the objective function for these optimizations is (50) Here, N B = ω B /h B (i.e., N B + 1 is the number of sampling points), ω k = kh B , and h B is the sampling interval. Except for the different choice of the objective function, the rest of the optimization procedure is the same as that described in Sec. V. In optimization runs that minimized J B , we considered four bandwidth values ω B /2π = {25, 50, 75, 100} MHz and used the fixed sampling interval h B /2π = 1 MHz.
For illustration purposes, we use a logarithmic measure of average squeezing, Q − (ω B ) = 10 log 10 P − (ω B ), however, note that Q − (ω B ) = Q − (ω) . Table VII shows the best values of Q − (ω B ) for ω B /2π = 100 MHz, for both the CQFC network of two coupled OPOs and the single OPO network, obtained in optimizations with T u = 0.9, L = L in = 0.01, and various values of L tl = L out and x u . We see that the CQFC network of two coupled OPOs significantly outperforms the single OPO network in terms of the average squeezing generated over the 100 MHz bandwidth, especially for lower values of transmission line losses.  It is also interesting to examine the squeezing spectrum Q − (ω) generated under the optimal operation of either network when we minimize J B = P − (ω B ). Figure 17 shows Q − (ω) for both networks for various values of ω B , L tl = L out , and x u . Similarly to the results shown in Sec. VI (cf. Fig. 10), we find that the optimally operated single OPO network generates exactly the same Lorentzian squeezing spectrum for any choice of ω B . In contrast, the CQFC network of two coupled OPOs is capable of adapting the generated squeezing spectrum depending on the selected value of ω B and overall produces much higher squeezing bandwidth.

X. CONCLUSIONS
We modeled the squeezing spectrum of the output field of the CQFC network of two coupled OPOs and used a suite of global optimization methods to examine the limits to which this spectrum can be varied under conditions typical for tabletop experiments. We found that, in contrast to a single OPO, the CQFC network can utilize the interference between the fields in the plant OPO and the controller OPO to significantly modify the squeezing spectrum of the output field in response to the selected optimization objective. In particular, when the objective is to maximize the degree of squeezing at a highfrequency sideband ω opt , the CQFC network can operate in an optimal regime characterized by a high degree of cooperativity between the plant OPO and the controller OPO, as quantified by the flows of light between them and the correlation between the phase shifts φ 1 and φ 2 . In this operation regime, the optimized squeezing spectrum Q − (ω) of the CQFC network of two coupled OPOs has the minimum at ω = ω opt , while the minimum of the optimized spectrum of the single OPO network is always at zero sideband frequency.
For both types of optimization objectives considered in this work (maximizing the degree of squeezing at a selected sideband frequency and maximizing the average degree of squeezing over a selected bandwidth), the CQFC network of two coupled OPOs significantly outperforms a single OPO in terms of squeezing achieved under similar conditions, even with higher losses in the CQFC network due to additional components and transmission lines. Also, the CQFC network is more effective in terms of converting a higher pump power into a stronger squeezing. While this superior performance of the CQFC network of two coupled OPOs relies on a phasesensitive interference between multiple fields, we discovered, perhaps surprisingly, that squeezing generated by the optimally operated CQFC network is rather robust to small variations of phase parameters. This robustness can be attributed to the tendency of global optimization algorithms to avoid solutions that are overly sensitive to small parameter variations, but the fact that such robust network configurations do actually exist is quite remarkable.
Overall, our results strongly indicate that CQFC networks provide a very effective tool for engineering quantum optical systems with new properties and unprecedented levels of performance. This work also demonstrates the usefulness of advanced optimization methods for analyzing and improving the performance of such networks.