Brian Legros
email: brian.legros@cnrs-orleans.fr

Luc Léger

Amer Chpoun

Mohamed Sellam

Azeddine Kourta
email: azeddine.kourta@univ-orleans.fr

Potential benefits of radial secondary injection of helium in dual-bell nozzles

Radial secondary injection of air in the extension section of dual-bell nozzles (DBN) have been the subject of recent studies. The technique showed promise in solving the early transition problem in DBNs, but also in reducing the side-loads generated during the transition phases. Because helium is used to pressurise rocket launcher's fuel tanks, the idea of using radial secondary injection of helium in the extension section of a DBN became of interest. This paper compares the impact of a radial secondary injection of helium in a DBN to a parametrical investigation performed with secondary injection of air. It was found that the use of helium resulted in a transition nozzle pressure ratio (NPR) and a retransition NPR increase of 5.5% and 11.4% respectively compared to test cases using secondary injection of air with a similar secondary pressure ratio. The lateral forces generated by the nozzle during the transition and the retransition phases were brought below 1% of the nozzle thrust by using helium or air as injectant gas. Moreover, the use of helium suppressed the hysteresis effect by causing an identical specific impulse curve trajectory during the transition phases and the retransition phases. Finally, the thrust jumps during transitions and retransitions were dramatically decreased, particularly with the use of helium, suggesting the promising potential of DBN on rocket launchers.

INTRODUCTION

The ever-increasing number of satellites to be launched into space have triggered numerous companies to emerge in order to get a market share. The rising number of competitors have pushed companies and state agencies to decrease their launching cost and to increase their launchers' efficiency to remain competitive. Launchers such as Ariane 5 use parallel staging which requires the first stage engine to be ignited before the lift off for reliability reasons [START_REF] Stark | Ariane 5 performance optimization using dual-bell nozzle extension[END_REF]. Therefore, the first stage engine is doomed to suffer from off-design operating conditions during the launcher's ascent, which induces high losses. Since the nozzle geometry is designed for a specified performance and the ambient pressure decreases during the launcher's ascent, the nozzle experiences different flow conditions: overexpansion, adapted flow and underexpansion (see Fig. 1). To prevent flow separation inside rocket nozzles at low altitudes (see Fig. 1) and therefore avoid flow instabilities and unsteady side-loads, the nozzle exit area is reduced, even though it implies lower performances at high altitudes. Ideally, the perfect nozzle should be able to increase its exit area during the launcher's ascent in order to adapt to the ambient pressure variations. Nevertheless, this solution would dramatically increase the complexity of the nozzle and decrease its reliability. Many nozzle concepts with altitude adaptive capabilities have been investigated in the past, yet this paper deals solely with the dual-bell nozzle concept. The dual-bell nozzle (DBN) concept was first mentioned by Swan [START_REF] Swan | The influence of design on the flight performance of rocket vehicles, with an analysis of the results of jet separation[END_REF] and consists of a converging/diverging nozzle, where the diverging section is com-posed by two successive nozzle profiles of different expansion ratios. The first part of the divergent, called the base nozzle, provides better performances at low altitudes with little risks of side-loads. The second part of the divergent, called the extension nozzle, provides better performances at high altitudes. The DBN geometry used in the present paper is illustrated in Fig. 2. The base nozzle and the extension nozzle are linked together through an inflection region. In sea-level mode, the flow detaches from the nozzle wall at the inflection point, providing a controlled and symmetric separation, therefore limiting the generation of lateral forces. As the launcher ascends in the atmosphere and the ambient pressure decreases, the flow will suddenly reattach the nozzle extension wall, moving the separation location at the nozzle exit. The attached flow in the entirety of the nozzle will provide a higher nozzle pressure ratio and therefore better performances in high altitudes. Fig. 3 shows schlieren images of the DBN in low-altitude mode and high-altitude mode during an experiment performed in the wind tunnel described in Sec. 2. In low-altitude mode, due to the high ambient pressure, the flow separates at the inflection point, which is the onset of the separation shock location [START_REF] Stark | Ariane 5 performance optimization using dual-bell nozzle extension[END_REF]. The separation shock propagates downstream until being reflected on the nozzle symmetry axis. The reflected shock (3) is then reflected on the exhaust plume boundary (2) further downstream. This process continues to form a succession of compression and expansion cells. In high-altitude mode, the flow is attached to the extension nozzle wall and the DBN is said to be full flowing. In this case, an expansion fan ( 5) is present at the inflection point. In Fig. 3, because the DBN extension profile geometry is characterised by a constant wall pressure, an internal recompression shock (6) emerges from the inflection point and directly interacts with the nozzle lip shock [START_REF] Ferrero | Dual-bell nozzle with fluidic control of transition for space launchers[END_REF]. The nozzle lip shock intersects with a Mach disc (8) at the triple point further downstream (see Fig. 3). The reflected shock (3) issued from the triple point is then reflected on the exhaust jet boundary [START_REF] Swan | The influence of design on the flight performance of rocket vehicles, with an analysis of the results of jet separation[END_REF]. This process also continues downstream, forming a succession of compression and expansion cells. The promising potential of dual-bell nozzles to reduce the cost of access to space fuelled an increasing number of investigations in many research centers [START_REF] Manuel | Critical assessment of dual-bell nozzles[END_REF][START_REF] Hagemann | Experimental and analytical design verification of the dual-bell concept[END_REF][START_REF] Kusaka | Experimental study on extendible and dual-bell nozzles under high altitude conditions[END_REF][START_REF] Dumnov | Dual bell nozzles for rocket engines of launch vehicle upper stages and orbital transfer vehicles[END_REF]. Optimisation codes developed to quantify the financial gain on an Ariane-5-like launcher mounted with a dualbell nozzle revealed a significant payload mass gain into geostationary transfer orbit of 490 kg [START_REF] Stark | Ariane 5 performance optimization using dual-bell nozzle extension[END_REF]. Recent analysis performed in [START_REF] Ferrero | Dual-bell nozzle with fluidic control of transition for space launchers[END_REF] for similar operating conditions indicated a payload mass gain of 1.5 tons using a DBN equipped with an optimised radial secondary injection in order to control the transition. The DBN appears to be a promising and easy-to-implement concept, with no cooling problems. However, three main issues have to be mindfully considered to make the DBN a viable concept: 1) Early transition and retransition, 2) Side forces generation, 3) Stability. Ideally, the transition from lowaltitude mode to high-altitude mode (and vice versa) must take place as close as possible to the optimum transition point (see Fig. 4) to provide the best performances. However, the natural transition and retransition in DBNs occur before the optimum transition point (as seen Fig. 4), causing the loss in specific impulse. Consequently, an efficient and reliable method has to be implemented to delay these critical phases. Besides, the shift from one mode to the other triggers dangerous side-loads caused by unsteady and asymmetrical flow separation in the extension section, which could be detrimental to the rocket integrity [START_REF] Nürnberger-Genin | Flow transition in dual bell nozzles[END_REF][START_REF] Génin | Side loads in subscale dual bell nozzles[END_REF]. The transition process in DBNs depends on several parameters that are still the subject of numerous studies: Reynolds Number and temperature, testing environment, nozzle geometry, ambient pressure fluctuations, etc. [START_REF] Génin | Side loads in subscale dual bell nozzles[END_REF][START_REF] Verma | Gas density effects on dual-bell transition behavior[END_REF][START_REF] Zmijanovic | Assessment of transition regimes in a dual-bell nozzle and possibility of active fluidic control[END_REF][START_REF] Génin | Experimental investigation of the inflection geometry on dual bell nozzle flow behavior[END_REF][START_REF] Génin | Experimental and numerical study of heat flux in dual bell nozzles[END_REF].

Separated flow region

Shocks

M>1

Shear layer Expansion waves

Adapted nozzle flow

Overexpanded nozzle flow

Underexpanded nozzle flow

Shear layer Expansion waves

The extension nozzle profile type in DBNs is among the most important parameters in regard to the transition behaviour [START_REF] Nasuti | Role of wall shape on the transition in axisymmetric dual-bell nozzles[END_REF]. There are three different types of extension profiles: Negative wall Pressure gradient (NP), Positive wall Pressure gradient (PP), and Constant wall Pressure (CP). It was shown in [START_REF] Tomita | Investigation on characteristics of conventional-nozzle-based altitude compensating nozzles by cold-flow tests[END_REF] that the separation location depends on the nozzle pressure ratio in NP extension profiles and was well predicted by flow separation criteria such as the Schmucker's criterion. However, the progressive movement of the separation front is the cause of important side-loads [START_REF] Hagemann | Experimental and analytical design verification of the dual-bell concept[END_REF][START_REF] Tomita | Investigation on characteristics of conventional-nozzle-based altitude compensating nozzles by cold-flow tests[END_REF]. In the case of PP extension profiles and CP extension profiles, a sudden transition is observed experimentally [START_REF] Hagemann | Experimental and analytical design verification of the dual-bell concept[END_REF][START_REF] Tomita | Investigation on characteristics of conventional-nozzle-based altitude compensating nozzles by cold-flow tests[END_REF] and numerically [START_REF] Nasuti | Numerical study of transition between the two operating modes of dual-bell nozzles[END_REF]. The sudden movement of the separation front in the extension reduces the risks of side-loads generation during the transition and the retransition phases [START_REF] Hagemann | Experimental and analytical design verification of the dual-bell concept[END_REF][START_REF] Tomita | Investigation on characteristics of conventional-nozzle-based altitude compensating nozzles by cold-flow tests[END_REF][START_REF] Horn | Dual-bell altitude compensating nozzles[END_REF]. The extension length also plays an important role in DBNs as a shorter extension would decrease the lever arm, and consequently reduce the side-loads moment. However, a decrease in extension length reduces the hysteresis effect (difference between the transition NPR and the retransition NPR) of the DBN and therefore its stability [START_REF] Génin | Side loads in subscale dual bell nozzles[END_REF]. Numerous studies showed that DBNs were sensitive to the testing environment, namely whether the tests were performed in altitude chambers or at ambient conditions. The variation of the feeding total pressure to increase (resp. decrease) the nozzle pressure ratio during experiments significantly influences the behaviour of DBNs. It changes the Reynolds number during the experiment and modifies the effective wall shape of the inflection region [START_REF] Verma | Gas density effects on dual-bell transition behavior[END_REF][START_REF] Nasuti | Role of wall shape on the transition in axisymmetric dual-bell nozzles[END_REF][START_REF] Verma | Reynolds number influence on dual-bell transition phenomena[END_REF][START_REF] Martelli | Numerical parametric analysis of dual-bell nozzle flows[END_REF][START_REF] Stark | Flow separation in rocket nozzles under high altitude condition[END_REF].

Even though many of the previous studies brought to light important information regarding DBNs behaviour during the ascent and the descent phases, few have focused on the control of transition and retransition in DBNs. Flow control in rocket nozzles has been the subject of previous studies, whether to prevent flow separation in overexpanded nozzles, to control thrust vector, or to decrease side-loads amplitude on a launching pad [START_REF] Zmijanovic | Experimental-numerical parametric investigation of a rocket nozzle secondary injection thrust vectoring[END_REF][START_REF] Stark | Experimental study on rocket nozzle side load reduction[END_REF]. Flow control in dual-bell nozzles has been mainly investigated by using film cooling and secondary injection normal to the mainstream direction, principally to control the transition.

Investigations on film cooling in the vicinity of the inflection region showed a decrease in transition NPR in [START_REF] Schneider | Active control of dual-bell nozzle operation mode transition by film cooling and mixture ratio variation[END_REF], but an increase of it for Proshchanka et al. in [START_REF] Proshchanka | Control of operation mode transition dual-bell nozzles with film cooling[END_REF]. A decrease in thrust jump (an example of thrust jump is shown in Fig. 4) during the transition phases, and lower lateral forces were also reported. In [START_REF] Martelli | Numerical analysis of film cooling in advanced rocket nozzles[END_REF], Martelli et al. showed an effective capability of lowering the wall temperature at the cost of long-lasting lateral forces capable of exciting a nozzle structural mode.

The use of radial secondary injection showed that outstanding gains could be found even for relatively small secondary mass flow rates. A radial secondary injection of air behind the inflection region significantly increased the transition NPR and the retransition NPR [START_REF] Léger | Controlled flow regime transition in a dual bell nozzle by secondary radial injection[END_REF]. The secondary fluid acts as an obstacle to the upstream flow, triggering a forced and symmetrical flow separation in the vicinity of the inflection region. Using the radial secondary injection, the flow will remain detached on a longer NPR range and bring one towards the optimum transition point [START_REF] Zmijanovic | Assessment of transition regimes in a dual-bell nozzle and possibility of active fluidic control[END_REF][START_REF] Léger | Controlled flow regime transition in a dual bell nozzle by secondary radial injection[END_REF]. Furthermore, Léger et al. [START_REF] Léger | Experimental investigation of forced flow regime transition in a dual bell nozzle by secondary fluidic injection[END_REF] showed that radial secondary injection impacts the nozzle stability and it reduces side-loads generation. It was also found that secondary injection reduced the flipflopping behaviour of the shock system during the retransition phases, where the flip-flopping phenomenon corresponds to a high amplitude and unsteady displacement of the separation point.

In this paper, the authors study the potential benefits of helium radial secondary injection in a dual-bell nozzle compared to the secondary injection of air.

EXPERIMENTAL SETUP

The experiments were performed in the EDITH depressurized wind tunnel (WT) (see Fig. 5), the same test facility as described in [START_REF] Léger | Experimental investigation of forced flow regime transition in a dual bell nozzle by secondary fluidic injection[END_REF]. EDITH is one of the rare nozzle testing facility in Europe offering nozzle exit pressure variations. Dry air from 320 litre tanks and pressurised at 300 bar is regulated to 3.5 bar to generate the DBN feeding total pressure in the stagnation chamber. The air from the stagnation chamber travels through the dual-bell nozzle and exits in the depressurised WT test section. The pressure inside the test section is controlled by a valve located in the WT diffuser, upstream of an overall 345 kW pumping group. The control of the position of the valve allows the decrease and the increase of pressure inside the WT test section. During the experiments, the dualbell nozzle is mounted on a force balance designed by the authors and which moves freely along the (x,y,z) orthogonal axis. The thrust and lateral forces are measured by the force balance using four HBM S2 strain-gauge force transducers, whose signals are amplified to a 0-10 V range before being acquired by the SCXI-1140 cards, at the rate of 1 kHz. The HBM S2 is a force sensor with an accuracy class of 0.05 which provides high accuracy measurements with a sensitivity of 2mv/V. Two transducers of 200 N located on both sides of the force balance were used to measure the vertical force component; one 200 N transducer measured the nozzle thrust; and a 20 N transducer measured the lateral forces. Previous studies for the validation of the force balance measurements showed a 0.24 N standard deviation for the force measurements on the longitudinal axis, 0.06 N on the vertical axis and 0.17 N on the lateral axis.

Dual-bell Nozzle

Wind tunnel The base nozzle was designed as a truncated ideal contour (TIC) for an exit Mach number of 3.2, and the extension nozzle was built as a constant pressure (CP) extension profile. To actively control the flow transition and retransition behaviour, a secondary flow of helium or air was radially injected in the extension part of the DBN. The 0.2 mm width injection slot was made possible by manufacturing the DBN in two parts, and the slot was placed 8 mm downstream of the inflection point. This technique also allowed the creation of a settling chamber for the secondary injection (called cavity, and visible in Fig. 6), to provide an axisymmetric, homogeneous secondary injection.

Schlieren visualization system

Secondary injection pipes

In this study, a test campaign has been carried out to compare the influence of air and helium as injectant gas on the DBN transition behaviour. In this case, several secondary injection pressures were used for air, ranging from 0.37 bar to 2.17 bar. Only one secondary injection pressure was used for helium: 0.31 bar. For all experiments, the DBN feeding total pressure was kept constant at 3.5 bar and the ambient pressure was varied using the main exit valve located upstream of the pumping group to trigger the transitions and the retransitions in the DBN. When secondary injection was used, the injection pressure was kept constant throughout the experiment. The ambient pressure and the secondary feeding total pressure inside the cavity were measured using Kulite XCQ-062 pressure transducers of 0-100 kPa range. For secondary injection pressures above 90kPa, the pressure transducer inside the cavity was replaced by a Kulite XCQ-093 of 0-350 kPa with a sensitivity of 28.743 mV/bar. Geometrical details of the DBN profile obtained by an in-house code based on the method of characteristics are presented in Fig. 2, and the reader can find further details in [START_REF] Léger | Experimental investigation of forced flow regime transition in a dual bell nozzle by secondary fluidic injection[END_REF]28].

RESULTS AND DISCUSSIONS

Helium and air secondary injection effects on dual-bell nozzle behaviour

This section exposes the impact of air and helium as secondary injectant gas on a dual-bell nozzle transition behaviour. The study compares the results obtained using one secondary injection pressure for helium to several injection pressures for air. During each experiment, the feeding total pressure was kept constant and the ambient pressure in the wind tunnel was repetitively increased and decreased to trigger the transitions and the retransitions inside the DBN. For each experiment, the chosen secondary injection pressure was kept constant. Several secondary injection pressures were considered. The secondary pressure values (P i ), the secondary pressure ratios (SPR=P i /P 0 ), the theoretical secondary mass flow rates ( ṁi ) and the secondary mass flow rate ratios (ϕ ṁ = ṁi / ṁ) are detailed in Tab. 1. Here, ṁ corresponds to the mainstream, theoretical mass flow rate of 0.187 kg/s. Fig. 7 shows the evolution of the measured parameters during a typical experiment using a constant air SPR of 0.15. At SPR = 0.15 and below NPR = 17.76, the dual-bell nozzle is in the low-altitude mode, where the flow detaches from the wall in the vicinity of the inflection region. When the NPR (red curve in Fig. 7) is increased from its lowest values, the specific impulse generated by the nozzle increases until a sudden drop in specific impulse occurs at NPR = 17.77. At NPR = 17.77, the transition takes place and the flow reattaches the extension wall. During this phase, the DBN switches from the low-altitude mode, to the high-altitude mode. Reciprocally, decreasing the NPR from the high-altitude mode will eventually trigger the retransition process, switching the DBN working mode from high-altitude, to lowaltitude. Previous studies on the present DBN revealed important side-loads peaks during transition and retransition phases ( [START_REF] Zmijanovic | Assessment of transition regimes in a dual-bell nozzle and possibility of active fluidic control[END_REF][START_REF] Léger | Controlled flow regime transition in a dual bell nozzle by secondary radial injection[END_REF][START_REF] Léger | Experimental investigation of forced flow regime transition in a dual bell nozzle by secondary fluidic injection[END_REF]), however, Fig. 7 shows that side-loads generation can be mostly suppressed with the use of radial secondary injection, as the peaks are hardly identifiable.

Instantaneous data during a transition and a retransition for several SPR cases are displayed in Fig. 8 and Fig. 9. Both figures display, on top, the specific impulse as a function of nozzle pressure ratio for the DBN without injection slot (called smooth DBN and visible in Fig. 6), the DBN with two secondary injection of air and one secondary injection of helium. On the bottom, Fig. 8 and Fig. 9 show the side-loads-to-thrust ratio. The instanta-neous data measurements revealed that, not only the use of secondary injection of air delays the transition NPR (NPR trans ) and the retransition NPR (NPR retrans ), but it also reduces the side-loads generated during the transition from one mode to the other. As shown in Fig. 8, for SPR = 0.15, the side-loads peaks could not be identified as the signal falls within the measurements noise. The use of secondary injection of helium significantly delays the NPR trans and the side-loads peaks fall within the measurement noise. During the descent phases, the use of secondary injection of air and helium increases the NPR retrans and outstandingly reduces the lateral forces generated (see Fig. 9). The evidence of the previously described behaviour is investigated in more details below with averaged data measurements with their standard deviation.

Fig. 10 shows the evolution of the NPR trans and NPR retrans as a function of SPR. The configuration associated with a SPR of 0 on Fig. 10, Fig. 12 and Fig. 13, corresponds to the case with DBN with secondary injection slot but without any injection. The experiments performed with radial secondary injection in the extension profile confirm that the NPR trans and the NPR retrans are greatly influenced by the SPR and the properties of the injectant gas (see Fig. 10). The use of air as injectant gas revealed an increase in NPR trans when the SPR was increased. The maximum NPR trans reached 18.39 and was achieved at SPR = 0.25, corresponding to an increase of 15.5% compared to the SPR = 0 case and an increase of 23.9% compared to the smooth nozzle configuration. Further increase in SPR showed a decrease in NPR trans . The largest decrease in NPR trans following the maximum value achieved was reached at SPR = 0.62 and amounted to 5.7%, bringing down the NPR trans to 17.34. A similar behaviour was identified for the retransition process, as the NPR retrans increased with an increase in SPR up to SPR = 0.36. The maximum NPR retrans reached 17.45, corresponding to an increase of 17.8% and 20.1% compared to the SPR = 0 case and smooth case, respectively. Similarly, further increase in SPR showed a decrease in NPR retrans . The largest decrease in NPR retrans following the maximum value achieved was reached at SPR = 0.62 and amounted to 0.9%, bringing down the NPR retrans to 17.29. Fig. 10 shows that delaying the NPR trans and the NPR retrans is limited as there exists a SPR from which an increase in secondary injection pressure causes a decrease in NPR trans and NPR retrans . High SPR might change the physics involved in the flow and further investigations using CFD should provide explanations to the previous phenomenon. In the helium configuration, only one secondary pressure ratio was considered: SPR = 0.09. It was found that the NPR trans and the NPR retrans were further increased compared to the results obtained with air for a similar SPR. Compared to the air SPR = 0.11 configuration, the experiments per- However, in this case, the NPR trans was less than the NPR retrans . The previous remark does not imply that the hysteresis effect is reversed, it indicates that the hysteresis effect is suppressed. Fig. 11 shows the specific impulse as a function of NPR for a configuration using air and a configuration using helium during a transition and a retransition phase. It shows that the specific impulse trajectory during the transition phase is relatively close to the specific impulse trajectory during the retransition phase in the helium configuration compared the configuration using air. Because the NPR trans (resp. NPR retrans ) is taken when the specific impulse jump takes place during increasing (resp. descending) NPR phases, the superimposed curves imply higher NPR retrans than NPR trans .

It is worth mentioning that the secondary injection of helium and air for SPR = 0.11 were not made at sonic speed due to the low feeding injection pressure. Further experiments using higher SPR for the helium configuration might lead to a better understanding of the aforementioned phenomenon.

Regarding the air configuration, the hysteresis was the largest for the SPR = 0 configuration, reaching 7%, then decreased to 2.7% for SPR = 0.11, and progressively increased along with the secondary injection pressure to hit a maximum of 5.7% for SPR = 0.2. A further increase in SPR induced a decrease in stability as the hysteresis became smaller, dropping to 0.3% for the highest SPR. The influence of the subsonic jet at SPR = 0.11 may have played a role in the early retransition in this specific case, as the cavity is free to suck or blow air freely in the vicinity of the injection slot region. Fig. 12 shows that the side forces were reduced by using radial secondary injection in both configuration: air and helium. At SPR = 0.09 for helium and for SPR = 0.15 and above for air, the side-loads generated during the transitions and the retransitions fall within the balance measurement noise (± 0.5 N) and the force balance could not measure any side-load peaks. At SPR = 0.09 for helium and at SPR ≥ 0.11 for air, the side-loads during transition dropped below 1% of the nozzle thrust, against 2.4% of the nozzle thrust in the smooth DBN case. Secondary injection of air brought side-loads during retransition below 1% of the nozzle thrust for SPR = 0.15 and above, against 3.4% in the smooth DBN configuration, providing over three times lower lateral forces magnitude compared to the smooth nozzle configuration. Secondary injection of helium at SPR = 0.09 also brought side-loads generation during retransition phases below 1% of the nozzle thrust. From SPR = 0.36 and above, a slight increase in lateral forces is visible during both, transition and retransition phases. This rise in lateral forces also falls within the measurement noise.

The influence of secondary injection on thrust jump during the transitions and the retransitions is showed in 13. The test series revealed that the air SPR did not present any significant impact on the thrust jump magnitude during transition and retransition between SPR = 0.15 and SPR = 0.25 as it stays fairly constant, though much less important than in the smooth nozzle configuration, which showed 2.2 times higher values. For SPR = 0.11, jumps magnitude during both, the ascent and the descent, are slightly higher than for the other cases. The latter observation might originate from the non-sonic secondary injection jet at SPR = 0.11 which can be a source of instability in the flow and a more thorough investigation is necessary to validate this hypothesis. Above SPR = 0.25, a further decrease in thrust jump is visible as the secondary injection pressure is increased, though the maximum decrease remains fairly close to the values between SPR = 0.15 and SPR = 0.25. The maximum thrust jump decrease during transition using secondary injection of air led to 3.2 times lower values than measured in the smooth DBN configuration. Fig. 13 also reveals that no significant differences were measured between the SPR = 0 configuration and the different injection cases during transition. Nonetheless, the presence of air secondary injection reduced the thrust jump during retransition up to 67.6% compared to the SPR = 0 case, and 73.8% (a factor of 3.8) compared to the smooth nozzle configuration. The effect of helium secondary injection during transition and retransition was even more significant as the decrease in thrust jump compared to the smooth nozzle configuration exceeded 86% during the ascent and the descent phases. The significant error bars visible on Fig. 13 emanates from the difficulty to clearly identify the thrust jump magnitude in this configuration, especially because of the force balance measurements limitations.

CONCLUSIONS AND SUMMARY

This paper focused on comparing the impact of air and helium radial secondary injection on the DBN transition behaviour. Several air secondary pressure ratios (SPR) ranging from 0.11 to 0.62 were considered and one SPR for helium was used: 0.09. Compared to the smooth nozzle configuration (nozzle without secondary injection slot), the transition NPR was delayed at maximum by nearly 23.9% using an air SPR of 0.25, against an increase of 22.5% for a helium SPR of 0.09. Compared to the smooth nozzle configuration, the retransition NPR was increased at maximum by 20.1% for an air SPR of 0.36, against 28.6% for a helium SPR of 0.09. The several experiments using different air SPR showed that the transition and retransition NPR delay were limited. The side-loads were decreased below 1% of the nozzle thrust during the transition and the retransition phases using a secondary pressure ratio of air and helium of only 0.15 and 0.09, respectively. Moreover, secondary injection of air (resp. helium) reduced the thrust jump during transition up to 69.2% (resp. 86.3%) and during retransition up to 73.8% (resp. 86.8%) to the smooth nozzle configuration. A complete interpretation of the present results necessitates further numerical work and a parametrical investigation of the helium secondary injection will be an interesting aspect of DBNs regime transition control. Nonetheless, the previous results would be extremely suitable should dual-bell nozzles with secondary injection be used in any current or next generation of launchers.

Figure 1 :

 1 Figure 1: Conventional bell nozzle working modes.

Figure 2 :

 2 Figure 2: Geometrical parameters of the DBN profile used during the test campaign.

Figure 3 :

 3 Figure 3: DBN in low-altitude mode (NPR = 15) on the top and high-altitude mode (NPR = 20) on the bottom.

Figure 4 :

 4 Figure 4: DBN efficiency compared to conventional nozzles.

Figure 5 :

 5 Figure 5: EDITH wind tunnel testing facility.

Figure 6 :

 6 Figure 6: A) DBN without secondary injection (smooth DBN), B) DBN used for experiments with secondary injection, C) Exploded view of the DBN used for experiments with secondary injection, D) DBN mounted with secondary injection pipes in the EDITH wind tunnel.

Figure 7 :

 7 Figure 7: Evolution of the specific impulse, nozzle pressure ratio (NPR), side-loads, and secondary pressure ratio (SPR) as a function of time during an experiment of DBN with an air SPR = 0.15.

Figure 8 :

 8 Figure 8: Instantaneous measurements of specific impulse (on top) and side-loads-to-thrust ratio (on bottom) during one transition phase for the smooth DBN, DBN with secondary injection of air and DBN with secondary injection of helium.

Figure 9 :

 9 Figure 9: Instantaneous measurements of specific impulse (on top) and side-loads-to-thrust ratio (on bottom) during one retransition phase for the smooth DBN, DBN with secondary injection of air and DBN with secondary injection of helium.

Figure 10 :

 10 Figure 10: Transition and retransition NPR as a function of secondary pressure ratio. Error bars show the standard calculated for each configuration.

Figure 11 :

 11 Figure 11: Specific impulse during a transition phase and a retransition phase as a function of NPR for secondary injection of air at SPR = 0.11 (top) and helium at SPR = 0.09 (bottom).

Figure 12 :

 12 Figure 12: Side-loads-to-thrust ratio measured during transitions and retransitions as a function of secondary pressure ratio. Error bars show the standard deviation calculated for each configuration.

Figure 13 :

 13 Figure 13: Thrust jump divided by the thrust during transitions and retransitions as a function of secondary pressure ratio. Error bars show the standard deviation calculated for each configuration.

Table 1 :

 1 Secondary injection pressure, pressure ratio, mass flow, and mass flow ratio used in the test campaign. revealed an increase in NPR trans and NPR retrans of 5.5% and 11.4% respectively. The NPR trans and the NPR retrans reached 18.19 and 18.68 respectively.

	P i [bar] 0.37	0.51	0.69 0.88	1.26	1.63	2.17	0.31 (Helium)
	SPR	0.11	0.15	0.2	0.25	0.36	0.47	0.62	0.09
	ṁi [g/s] 1.98	2.73	3.7	4.72	6.73	8.75	11.63 4.79
	ϕ ṁ	0.011 0.015 0.02 0.025 0.036 0.047 0.062 0.026
	formed with helium

ACKNOWLEDGMENT

We would like to gratefully acknowledge the laboratory of Excellence CAPRYSSES framework and the financial support of this study from Grant No. ANR-11-LABX-0006-01 of the Investissements d'Avenir LabEx CAPRYSSES. would also like to thank Nicolas Gouillon for his technical assistance.


Mathilde Fruchart
email: mathilde.fruchart@univ-lille.fr.

Sylvia Pelayo

Christian Vilhelm

Antoine Lamer

Jupyter Notebooks for Introducing Data Science to Novice Users

Keywords: Data Science, Education, Social media, Programming 1

Data science is a bridge discipline involving computer science, statistics, and knowledge of the health field. We developed a Jupyter Notebook to enable novice users to easily and autonomously analyze data from social networks. We conducted an experimentation with non-programmer students. They had to adapt a R Notebook and complete 14 questions and to perform descriptive analyses. The average rate of correct answers was 90.7. Jupyter Notebook enabled novice users to easily and autonomously analyze data from Twitter.

Introduction

Data science may seem to be restricted to expert users because it requires the mastery of many technologies and knowledges [START_REF] Parker | Ten simple rules for starting (and sustaining) an academic data science initiative[END_REF]. Analyses are generally programmed from scripts consisting of a series of instructions, accompanied by comments, and developed with an integrated development environment (IDE) or in command line. Notebooks are another approach to programming and sharing results and were designed to improve the reproducibility of the studies [START_REF] Kluyver | Jupyter Notebooks -a publishing format for reproducible computational workflows[END_REF]. The final notebook provides a full and comprehensive analysis report. We developed a R notebook to enable novice users to easily and autonomously analyze data from social networks. We assessed their ability to use it and captured their feedback on the value of such a tool for their daily practice.

Methods

Jupyter Notebook is an interactive computational laboratory notebook, which can work with code in many different programming languages such as Python, Java, R, or Julia [START_REF]Project Jupyter[END_REF]. Jupyter Notebook allows for the smooth integration of code and narrative text (in Markdown syntax) into a single document that can be executed and edited on the fly. We experimented the Jupyter Notebook with students from two master's degree formations not specialized in Data Science and without programming skills.

Wedistributed a R Notebook with datasets from Twitter. had to execute the Notebook and to adapt the script in order to obtain additional results.

We assessed the ability of the students to use the Jupyter Notebook in checking if they were able to perform the 14 tasks [START_REF] Fruchart | Twitter analysis for novice users, GitLab[END_REF]: (i) instructions in relation to the notebook (modify and complete cells, execute cells, create new cells, export in pdf); (ii) R code 1 st level (from an existing code, compute a new variable, change a color in a graphic, change a variable in a graphic); (iii) R code 2 st level (find a instruction already implemented and modify it to produce a new result).

Results

We have delivered 20 health topics to 44 students. For each topic, we extracted 1000 tweets. Average (SD) rate of correct answers was 90.7 (13.3) for all the questions. Questions with the poorest scores were about (i) the modification of the required dataframe and the variable for the generation of a barplot, (ii) the updating of an instruction to display mentions instead of hastags and the adaptation of the title of the graphic, (iii) the modification of cells in Markdown format or (iv) the creation of two new cells for code or Mardown. The main difficulties that emerged were ensuring that users executed all the code cells, and that they had set up their cells in Markdown or code according to their needs. On the contrary, we did not report any major error on syntax problems which we usually encounter when learning a new language.

Discussion and conclusion

We developed a R notebook to enable novice users to easily and autonomously analyze data from Twitter. We tried it out with a class of 44 students without programming skills. Compared to Rmarkdown, which also integrate code and narrative text, Jupyter Notebook prevents the need to work locally and to have to install the necessary libraries himself [START_REF] Baumer | R Markdown: Integrating A Reproducible Analysis Tool into Introductory Statistics[END_REF]. IDE are rather designed for more advenced developers because they provide interesting features, such as code autocompletion, and displays the environment with all variables. However they could be more complex to handle for novice users than Jupyter notebook. Improvements would be to not display warnings and messages to provide training on manipulation of the notebook.


W C Smith

A Borgognoni

M Villiers

E Roverc'h

J Palomo

M R Delbecq

T Kontos

P Campagne-Ibarcq

B Douçot

Z Leghtas
email: zaki.leghtas@ens.fr

Spectral signature of high-order photon processes mediated by Cooper-pair pairing

Inducing interactions between individual photons is essential for applications in photonic quantum information processing and fundamental research on many-body photon states. A field that is well suited to combine strong interactions and low losses is microwave quantum optics with superconducting circuits. Photons are typically stored in an LC circuit, and interactions appear when the circuit is shunted by a Josephson tunnel junction. Importantly, the zero-point fluctuations of the superconducting phase across the junction control the strength and order of the induced interactions. Superconducting circuits have almost exclusively operated in the regime where phase fluctuations are smaller than unity, and two-photon interactions, known as the Kerr effect, dominate. In this experiment, we shunt a high-impedance LC oscillator by a dipole that only allows pairs of Cooper pairs to tunnel. Phase fluctuations, which are effectively doubled by this pairing, reach the value of 3.4. In this regime of extreme fluctuations, we observe transition frequencies that shift non-monotonically as we climb the anharmonic ladder. From this spectroscopic measurement, we extract two-, three-and four-photon interaction energies of comparable amplitude, and all exceeding the photon loss rate. This work explores a new regime of high-order photon interactions in microwave quantum optics, with applications ranging from multi-photon quantum logic to the study of highly correlated microwave radiation.

Photons do not interact with each other in free space. In the quantum optical domain, they are typically brought into interaction by coupling them to atoms [START_REF] Chang | Quantum nonlinear optics -photon by photon[END_REF]. Recent advances have realized two-and three-photon interactions mediated by a dense gas of Rydberg atoms, demonstrating photon dimers and trimers [START_REF] Liang | Observation of three-photon bound states in a quantum nonlinear medium[END_REF], and photonic vortices [START_REF] Drori | Quantum vortices of strongly interacting photons[END_REF]. Reaching processes of higher order would find applications in multi-photon quantum logic [START_REF] Adhikari | Nonlinear optics quantum computing with circuit QED[END_REF] and the study of many-body photon states [START_REF] Greentree | Quantum phase transitions of light[END_REF][START_REF] Ma | A dissipatively stabilized Mott insulator of photons[END_REF][START_REF] Carusotto | Photonic materials in circuit quantum electrodynamics[END_REF], but has remained out of reach since it requires inducing even stronger interactions between photons.

In the field of microwave quantum optics with superconducting circuits [START_REF] Wallraff | Strong coupling of a single photon to a superconducting qubit using circuit quantum electrodynamics[END_REF][START_REF] Hofheinz | Synthesizing arbitrary quantum states in a superconducting resonator[END_REF][START_REF] Lang | Correlations, indistinguishability and entanglement in Hong-Ou-Mandel experiments at microwave frequencies[END_REF], the nonlinearity of the Josephson junction is employed to mediate interactions between photons. These photons are typically stored in an LC circuit [START_REF] Devoret | Quantum fluctuations in electrical circuits[END_REF] of angular frequency Ω = 1/ √ LC (L and C are the circuit inductance and capacitance, respectively). Superconductivity endows these circuits with low photon loss, and quality factors exceeding one million are routinely observed [START_REF] Place | New material platform for superconducting transmon qubits with coherence times exceeding 0.3 milliseconds[END_REF][START_REF] Ganjam | Surpassing millisecond coherence times in on-chip superconducting quantum memories by optimizing materials, processes, and circuit design[END_REF][START_REF] Kono | Mechanically induced correlated errors on superconducting qubits with relaxation times exceeding 0.4 milliseconds[END_REF]. When such a circuit is shunted by a Josephson junction, interactions between photons appear (Fig. 1). The Hamiltonian takes the form

Ĥideal = ℏΩâ † â -E J cos( φ -ϕ ext ) φ = η(â + â † ) , (1)
where ℏ is the reduced Planck constant and â is the photon annihilation operator. The interaction energy stems from the Josephson cosine potential with Josephson energy E J , and φ is the phase drop across the junction with η its zero-point fluctuations. The loop formed by the oscillator inductance and the junction is threaded by external magnetic flux denoted ϕ ext . Expanding the cosine into its Taylor series reveals the various interaction processes. For example, the η(â + â † ) 4 term yields a two-photon interaction term â † 2 â2 corresponding to the Kerr effect. A celebrated success of microwave quantum optics was the first realization of a Kerr interaction that exceeded the photon loss rate, demonstrating the collapse of a coherent state into multi-component Schrödinger cat states [START_REF] Kirchmair | Observation of quantum state collapse and revival due to the single-photon Kerr effect[END_REF]. In this letter, we address the problem of inducing higher-order processes of the form â † n ân where n = 3, 4, and beyond.

The relative strength of multi-photon processes is governed by the dimensionless quantity η. It can be expressed as η = πZ/R Q , where Z = L/C is the LC circuit impedance, and R Q ≈ 6.4 kΩ is the superconducting resistance quantum [START_REF] Girvin | Circuit QED: Superconducting qubits coupled to microwave photons[END_REF]. The n-photon process â † n ân has a strength that scales as (η n /n!) 2 , and therefore, its activation for large n requires entering the regime η > 1, or equivalently Z ≫ R Q . However, fabricating an LC oscillator with a characteristic impedance exceeding the superconducting resistance quantum is challenging, although progress is being made with the development of superinductances [START_REF] Manucharyan | Fluxonium: Single Cooper-pair circuit free of charge offsets[END_REF][START_REF] Peruzzo | Surpassing the resistance quantum with a geometric superinductor[END_REF][START_REF] Pechenezhskiy | The superconducting quasicharge qubit[END_REF][START_REF] Rieger | Granular aluminium nanojunction fluxonium qubit[END_REF]. Recently, η ≈ 1 was achieved in a planar coil resonator, and emissions of k-photon bunches (k = 1 to 6) were observed by activating the process âk + â † k with a voltage-biased junction [START_REF] Ménard | Emission of photon multiplets by a dc-biased superconducting circuit[END_REF].

Another route towards large phase fluctuations is to replace the Josephson junction, that allows Cooperpair tunneling, by a dipole that only allows pairs of 1) with parameters in Tab. I] and its linear equivalent (dashed lines) as a function of the superconducting phase difference at ϕext = π. Crucially, the adjacent transition frequency shift δn from the bare frequency ω0 alternates in sign when ascending the ladder. (c) Photon process diagrams for the first four interaction orders and corresponding interaction energies Jn.

Cooper pairs to tunnel [Fig. 1(a)]. In the basis of tunneled Cooper-pair number N , the tunneling operator is transformed as

1 2 N (|N ⟩ ⟨N + 1| + |N + 1⟩ ⟨N |) → 1 2 N (|N ⟩ ⟨N + 2| + |N + 2⟩ ⟨N |).
Equivalently, in the conjugate phase representation φ, cos φ → cos 2 φ [START_REF] Douçot | Pairing of Cooper pairs in a fully frustrated Josephson-junction chain[END_REF][START_REF] Ioffe | man, Possible realization of an ideal quantum computer in Josephson junction array[END_REF]. Denoting φ = 2 × φ, we see that phase fluctuations are effectively doubled: η = 2 × πZ/R Q [START_REF] Smith | Superconducting circuit protected by two-Cooperpair tunneling[END_REF]. In the extreme regime of η > 1, we further require that E J ≪ ℏΩ [Fig. 1(b)], so that the Josephson cosine potential primarily induces n-photon interaction processes â † n ân [Fig. 1(c)].

In this experiment, we implement a superconducting LC oscillator shunted by a two-Cooper-pair tunnel element. We place ourselves in the unexplored regime where the tunneling energy is smaller than the oscillator transition energies, and the zero-point phase fluctuations exceed one: E J /ℏΩ = 0.28 and η = 3.4. We measure the first four transition energies of our device, and find that unlike a Kerr resonator, they do not follow a monotonic trend. Instead, we observe an alternation of the sign of the oscillator frequency shift for each added photon. From this spectroscopic signature, we extract two-, three-, and four-photon interaction processes of amplitudes greater than 70 MHz, that alternate in sign, and far exceed the transition linewidths of 200 kHz. The realization of these unprecedentedly high-order photon pro-FIG. 2. Signatures of higher-order photon processes. (a)-(b) Simulated Wigner quasiprobability distribution representing the initial time evolution of the coherent state |α⟩ with α = 1.7 for small and large quantum phase fluctuations. This value of α was chosen to match η/2 for η = 3.4 (see text). As quantum phase fluctuations increase, the evolution remarkably transitions from diffusive to nondiffusive. Simulated transition frequency shifts δn = ωn -ω0, where ωn is the transition frequency between energy levels n and n + 1, are represented versus photon number at the starred external flux value (c)-(d), and versus external flux (e)-(f). We observe the transition from an ordered to an alternating arrangement that asymptotically approaches a Bessel function (dashed lines). Simulations correspond to numerical diagonalization and time propagation of Eq. (1) with parameters in Tab. I. cesses opens avenues in microwave quantum optics, such as multi-photon quantum logic [START_REF] Adhikari | Nonlinear optics quantum computing with circuit QED[END_REF], the study of manybody photon states [START_REF] Ma | A dissipatively stabilized Mott insulator of photons[END_REF], or the processing of protected qubits [START_REF] Cohen | Degeneracy-preserving quantum nondemolition measurement of parity-type observables for cat qubits[END_REF].

We proceed to the analysis of the ideal Hamiltonian in Eq. (1) in the unexplored regime E J ≪ ℏΩ and η > 1 (Fig. 1), and for simplicity, we set ϕ ext = 0. Note that the cosine in Eq. (1) may be decomposed as cos η(â

+ â † ) = 1 2
Dη + D-η , where Dη = exp iη(â + â † ) is the displacement evolution operator. Remarkably, this evolution operator that usually results from the integration of a linear Hamiltonian ∝ â + â † over time enters the Hamiltonian directly. As a consequence, even at short times, a quantum state evolving under this Hamiltonian in Eq. (1) will be displaced across phase-space by ±η. The effect is particularly striking when initializing the system in a coherent state of amplitude α = η/2 [START_REF] Cohen | Degeneracy-preserving quantum nondemolition measurement of parity-type observables for cat qubits[END_REF], so that |±α⟩, that are distant by 2α is phase-space, are directly coupled through cos η(â + â † ) [Fig. 2(b)]. This is qualitatively different from the familiar diffusive-like evolutions resulting from low-order photon interactions [START_REF] Kirchmair | Observation of quantum state collapse and revival due to the single-photon Kerr effect[END_REF] [Fig. 2(a)].

We now express Ĥideal in Eq. ( 1) in terms of n-photon interaction processes [Fig. 1(c)]. We start by expanding the cosine into a normal ordered Taylor expansion. Then, by virtue of the rotating wave approximation (RWA), we neglect non-particle number conserving terms. We arrive at [26] :

Ĥideal ≈ ℏω 0 â † â + n≥2 J n â † n ân , (2)
where the n-photon interaction energy takes the form

J n = -E J e -η 2 /2 (-1) n (η n /n!) 2 ,

and the renormalized frequency is ω

0 = Ω + J 1 /ℏ. Note that J n /J n-1 = -(η/n) 2 ,
and hence the interaction strength is maximal for the integer order closest to η. The eigenenergies of Eq. ( 2) are

E n = nℏω 0 + n k=2 n! (n-k)! J k ,
and the experimentally accessible quantities are the transition frequencies

ω n = (E n+1 -E n )/ℏ.
We introduce the transition frequency shift in the presence of n photons as δ n = ω n -ω 0 [Fig. 1(b)] and we find :

δ n = n+1 k=2 k n! (n + 1 -k)! J k /ℏ . (3)
In the familiar situation of the Kerr oscillator where . This resembles the oscillatory nonlinearity predicted in a resonator containing a phase-slip element [START_REF] Hriscu | Model of a proposed superconducting phase slip oscillator: A method for obtaining few-photon nonlinearities[END_REF].

η ≪ 1, J 2 ≈ -E J η 4 /
Our circuit implementation of Ĥideal in Eq. ( 1) is depicted in Figs. 3(a from the capacitance between the two wires linking the oscillator to the tunneling element, resulting in charging energy ϵ C /h = 5.65 GHz. Tunneling occurs through a socalled Kinetic Interference coTunneling Element (KITE) [START_REF] Smith | Superconducting circuit protected by two-Cooperpair tunneling[END_REF][START_REF] Bell | Protected Josephson rhombus chains[END_REF]. It consists of two parallel arms that form a loop threaded by an external flux θ ext . Each arm contains a small junction of Josephson energy

E ± J = E J (1 ± ϵ) and charging energy E ± C = E C /(1 ± ϵ)
, with E J /h = 3.92 GHz, E C /h = 6.23 GHz and fabrication uncertainty results in a small asymmetry factor ϵ = 0.032. Each small junction is placed in series with 6 large junctions of a total inductive energy ϵ L /h = 6.73 GHz. This nonvanishing series inductance induces the tunneling of pairs of Cooper-pairs. Finally, this circuit is inductively coupled to a resonator for readout [START_REF] Smith | Quantization of inductively shunted superconducting circuits[END_REF].

The circuit parameters quoted above are extracted by fitting a three-mode circuit model [26] (discarding the readout mode) to two-tone spectroscopy data at various flux biases (θ ext , φ ext ) [Figs. 3(c)-(f)]. This three-mode Hamiltonian is 2π periodic in (θ ext , φ ext ) and its spectrum possesses inversion symmetry about (θ ext , φ ext ) = (0, 0), (0, π), (π, 0), and (π, π) [START_REF] Smith | Magnifying quantum phase fluctuations with Cooper-pair pairing[END_REF]. Therefore, these four points are the vertices of a plaquette that constitutes the primitive cell of the circuit spectrum as a function of external flux. We acquire the circuit spectrum along the edges of this plaquette [diagrams in Figs. 3(c)-(f)]. At each bias point, we start by acquiring the reflection spectrum of the readout resonator [Figs. 3(g)-(j)]. We then set the readout tone on resonance, and sweep a probe tone over a broad spectral range [Figs. 3(c)-(f)]. When the probe hits a circuit transition from the ground state to the n-th excited state, the reflected readout signal is affected. We identify several transitions that are captured by a three-mode circuit model [26] (dashed lines). The feature near 4.4 GHz corresponds to a parasitic mode visible in electromagnetic simulations, and the one near 5.1 GHz corresponds to the readout resonator. Other non-captured features could arise from multi-photon transitions to higher excited states or junction array modes. The circuit parameters we extract from this fit are summarized in Tab. I.

Importantly-at sufficiently low frequencies (here be-low about 20 GHz)-our circuit is well described by a single-mode Hamiltonian that emulates the one in Eq.

(1) [26]. Discarding the readout mode and the two KITE self-resonant high-frequency modes, we focus on the LC oscillator of inductive energy E L and charging energy ϵ C . This oscillator is shunted by a KITE that renormalizes its frequency Ω and provides nonlinearity. In the regime ϵ L ≫ E J , the potential energy of one arm of the KITE traversed by a phase drop of φ takes the

form U ± (φ) ≈ -E ± J cos φ + E ± J 2
4ϵ L cos 2φ and higher harmonics have been neglected [26]. Biasing the circuit at θ ext = π, Cooper-pair tunneling across both arms interferes destructively, while cotunneling interferes constructively. Indeed, the potential energy of the KITE is

U + (φ) + U -(φ + π) ≈ -δE J cos φ + E J cos 2φ, with δE J = 2ϵE J and E J ≈ E 2 J /2ϵ L .
Conveniently, E J can be made smaller than the oscillator frequency Ω by adequately choosing E J and ϵ L . In summary, this yields an effective Hamiltonian for our circuit at θ ext = π of the form [26]:

Ĥcircuit = ℏΩâ † â + E J cos[2( φ -φ ext )] -δE J cos( φ -φ ext ) , (4)
where φ ext is the flux threading the loop formed by the KITE and the oscillator inductance. Note that Ĥcircuit corresponds to Ĥideal up to the perturbative term in δE J , with the correspondence φ = 2 φ and ϕ ext = 2φ ext + π.

Additionally, the J n obtained from Eq. ( 4) depend on external flux and contain an added contribution from the term in δE J .

We extract the parameters of the one-mode Hamiltonian in Eq. ( 4) by fitting this model to the measured transition energies at θ ext = π [Fig. 4(a)]. The resulting parameters are displayed in Tab. I. Two notable features are visible in the data. First, these transition frequencies vary in a ±500 MHz window-a ±17 % fraction of the central frequency Ω/2π = 2.86 GHz. This confirms that the flux dependent tunneling amplitude is a perturbation to the LC oscillator frequency, i.e. E J , δE J < ℏΩ. In particular, we find E J /h = 0.79 GHz and the perturbation δE J /h = 0.27 GHz. Second, the transition frequencies ω n between levels n and n + 1 are not ordered in n. Instead, they interlace as a function of φ ext , indicating that we have entered the regime of large phase fluctuations. For example, at φ ext = 0.2, ω 0 > ω 1 , ω 1 < ω 2 , and ω 2 > ω 3 [Fig. 4(b)]. From this measured spectrum, we compute the n-photon interaction strengths J n for n = 2, 3, and 4 [Fig. 4 In conclusion, this experiment explores a new regime of nonlinear microwave quantum optics where interactions between photons are so strong that second-, third-, and fourth-order processes are of comparable amplitude and largely exceed the photon decay rate. We access this regime with photons stored in a high impedance LC oscillator that is shunted by a two-Cooper-pair tunneling element, effectively boosting phase fluctuations. Two technical challenges must be met: the tunneling energy E J must be weaker than the oscillator energy ℏΩ, and the phase fluctuations η across the tunneling element must exceed one. We measure the first four transition frequencies of our circuit, and observe their interlacing versus flux. From these spectra, we extract E J /ℏΩ = 0.28 and η = 3.4. The observation of high-order photon processes could be extended in multiple directions. One is the study of the quantum dynamics and scattered radiation correlations of this system under the action of drives and dissipation [START_REF] Hriscu | Model of a proposed superconducting phase slip oscillator: A method for obtaining few-photon nonlinearities[END_REF]. Moreover, coupling our two-Cooperpair tunneling element to an array of resonators could induce high-order interactions between multiple modes, useful for the study of many-body photon states [START_REF] Greentree | Quantum phase transitions of light[END_REF][START_REF] Ma | A dissipatively stabilized Mott insulator of photons[END_REF]. Finally, applications are envisioned to process quantum information that is encoded non-locally over the phase space of an oscillator [START_REF] Cohen | Degeneracy-preserving quantum nondemolition measurement of parity-type observables for cat qubits[END_REF]. This section details the fabrication process we follow to produce the sample of this experiment.

Wafer preparation: The circuit is fabricated on a 430 µm-thick wafer of 0001-oriented, double-side epipolished Sapphire C. The sapphire wafer is initially cleaned through a stripping process in a reactive ion etching (RIE) machine, after which it is loaded into a sputtering system. After one night of pumping, we initiate an argon milling cleaning step, followed by the sputtering of 120 nm of niobium. Subsequently, we apply a protective layer of poly(methyl methacrylate) (PMMA A6), dice the wafer and clean the small chips in solvents. This is followed by a 2 min oxygen-stripping process and approximately 30 s of exposure to sulphur hexafluoride (SF6) in order to remove the oxide layer formed on the niobium during the stripping process.

Circuit patterning: We spin optical resist (S1805) and pattern the large features (control lines and readout resonator capacitor pads) using a laser writer. After development (MF319), we rinse in de-ionized water for 1 min, and etch the sample in SF6 with a 20 s over-etch. Finally, the sample is cleaned for 10 min in acetone at 50 • C.

Junction patterning: Next, we apply a bilayer of methacrylic acid/methyl methacrylate [MMA (8.5) MAA EL10] and poly(methyl methacrylate) (PMMA A6). The entire circuit (KITE, inductive shunt, and readout resonator) is patterned in a single e-beam lithography step. The development takes place in a 3:1 isopropyl alcohol (IPA)/water solution at 6 • C for 90 s, followed by 10 s in IPA. The undercut regions of the bilayer are cleaned by oxygen-stripping for 30 s.

Junction deposition:

The chip is then loaded in an e-beam evaporator. We start with a thorough argon ion milling for 2 min at ±30 • angles. We then evaporate 35 nm and 100 nm of aluminum, at ±30 • angles, separated by an oxidation step in 200 mbar of pure oxygen for 10 min.

Junction characteristics:

The Josephson junctions are all fabricated from Al/AlOx/Al in a single evaporation step, utilizing the Dolan bridge method. The e-beam base dose is set to 283 µC cm -2 , with an acceleration voltage of 20 kV and a lens aperture of 7.5 µm. Three types of junctions are fabricated. (i) Two small junctions are located in the KITE, with an area of 0.076 µm 2 . They are patterned with a dose factor of 0.9 and an undercut dose of 0.2, resulting in an inductance per junction of 42 nH. (ii) A total of 12 large array junctions are found within the KITE loop, along with a further 90 unshared shunting array junctions. These junctions have an area of 0.62 µm 2 , and are patterned with a dose factor of 0.9 and an undercut dose of 0.1, resulting in an inductance per junction of 4 nH. (iii) There are 29 larger array junctions that constitute the readout resonator inductance, 19 of which are shared. These junctions have an area of 3 µm 2 , and are patterned with a dose factor of 0.8 and an undercut dose of 0.1, resulting in an inductance per junction of 0.7 nH.

Sample mounting: The chip was subsequently glued with PMMA onto a PCB, wire-bonded and mounted into a sample holder. The device was then thermally anchored to the base plate of a Bluefors dilution refrigerator, surrounded by three concentric cans for magnetic and infrared shielding (outer: cryoperm, middle: aluminum, inner: copper). An optical micrograph of the circuit and SEM images of some junctions are shown in Fig. S5.

MATHEMATICAL DERIVATIONS

In this section we derive the expression of the n-photon interaction strength J n that enters Eq. ( 2), as well as Eq. ( 3) that relates J n to the transition frequency shift δ n .

We start from Hamiltonian (1) with ϕ ext = 0 for simplicity.

Note that the cosine may be decomposed as cos η(â + â † ) = 1 2 exp iη(â + â † ) + exp -iη(â + â † ) . Since the commutator [â † , â] commutes with both â and â † , we may use Glauber's formula as follows exp iη(â + â † ) = e -η 2 /2 e iηâ † e iηâ . We then expand each exponential into its Taylor series. Since we place ourselves in the regime ℏΩ ≫ E J , we perform the rotating wave approximation (RWA) by neglecting terms that do not conserve particle number, yielding cos φ

RWA = e -η 2 /2 +∞ n=0 (-1) n η n n! 2 â † n ân .
Inserting this last equation in Eq. ( 1) yields Hamiltonian (2) with :

J n = -E J e -η 2 /2 (-1) n (η n /n!) 2 .
We now turn to the derivation of Eq. ( 3). The eigenvectors of Hamiltonian ( 2

E n = nℏω 0 + n k=2 n! (n -k)! J k .
We define the transition frequency ω n = (E n+1 -E n )/ℏ, and δ n = ω n -ω 0 . Note that

E n+1 = (n + 1)ℏω 0 + (n + 1)!J n+1 + n k=2 (n + 1)! (n + 1 -k)! J k , so that ℏδ n = (n + 1)!J n+1 + n k=2 (n + 1)! (n + 1 -k)! - n! (n -k)! J k = (n + 1)!J n+1 + n! n k=2 k n! (n + 1 -k)! J k = n+1 k=2 k n! (n + 1 -k)! J k .

HARMONICS FROM SINGLE ARM KITE

In this section we demonstrate that a Josephson junction in series with a small inductance generates Josephson harmonics [START_REF] Willsch | Observation of Josephson harmonics in tunnel junctions[END_REF]. We carry out this simple derivation at zero frequency where capacitors may be disregarded. We consider the circuit depicted in Fig. S6, in the regime where the inductive energy ϵ L largely exceeds the Josephson energy E J . We denote ν = E J /ϵ L (here ν ≪ 1), φ the phase drop across the entire circuit, and φ J the phase drop across the junction alone. The potential energy U of this circuit is the sum of the inductive energy U L = 1 2 ϵ L (φ -φ J ) 2 and the Josephson energy U J = -E J cos φ J . This FIG. S6. Electrical circuit diagrams depicting the emergence of Josephson harmonics when a linear inductance is connected in series with a Josephson junction without its selfcapacitance (left). Imposing Kirchhoff's current law on the central node yields an effective potential containing both oneand two-Cooper-pair tunneling terms (right).

yields U = ϵ L 1 2 (φ -φ J ) 2 -ν cos φ J .
Also, φ J is tied to φ though Kirchhoff's current law E J sin φ J = ϵ L (φ -φ J ). This equation may be recast as φ J = φ -ν sin φ J , and to first order in ν, simplified to φ J = φ -ν sin φ + O(ν 2 ). Injecting this expression in the inductive energy yields

U L = 1 2 ϵ L ν 2 sin 2 φ + O(ν 3 ) and U J = -ϵ L ν cos φ + ν 2 sin 2 φ + O(ν 3 ) . Summing these two energies results in U = -ϵ L ν cos φ + 1 2 ν 2 sin 2 φ + O(ν 3 ) .
Finally, disregarding an irrelevant constant, we recover

U (φ) = -E J cos φ + 1 2 E 2 J 2ϵ L cos 2φ + ϵ L O (E J /ϵ L ) 3 .
(S5) The first term in this expression is the usual Josephson energy that allows Cooper-pair tunneling. The second term is the second Josephson harmonic, and is responsible for the tunneling of pairs of Cooper pairs (Fig. S6). The remaining perturbative term encompasses higher order harmonics.

MODEL REDUCTION

Outline of analysis

In this section, we describe the mathematical procedure for reducing the three-mode Hamiltonian of an LC oscillator shunted by a KITE to an effective one-mode Hamiltonian. The procedure that we use is not mathematically rigorous, yet it provides a reduced model that agrees both qualitatively and, to a lesser extent, quantitatively with the full circuit model (see Fig. S8). Improving the agreement and establishing a firm mathematical foundation will be the subject of future work [START_REF] Teufel | Adiabatic perturbation theory in quantum dynamics[END_REF]. We neglect the readout resonator and consider the circuit in Fig. S7. To realize the Hamiltonian in Eq. ( 4) at θ ext = π, we first place ourselves in the parameter regime E L ≪ ϵ C ∼ E C , which guarantees large phase fluctuations in the oscillator. Second, we require E J ≲ ϵ L ≪ E J /ε to lower the two-Cooper-pair tunneling energy below the bare Josephson energy E J while keeping single-Cooper-pair tunneling a perturbation. Finally, making the oscillator frequency the largest energy scale in the effective Hamiltonian (but still much smaller than the junction plasma frequency), we arrive at the parameter regime

E L ≪ ϵ C ∼ E C , E J ≲ ϵ L and ε ≪ 1.

FIG. S7. Electrical circuit diagram for the KITE (green)

shunted by a superinductance and a capacitance (blue). As in the main text, we have E ± J = (1±ε)EJ and

E ± C = EC /(1±ε).
Our analysis is organized as follows. We first write down the exact Hamiltonian for the circuit in Fig. S7 and move to a basis where the kinetic energy is diagonal and the modes are weakly coupled. Then, we derive the low-frequency component of the Hamiltonian evaluated at the equilibrium positions of the high-frequency modes. Finally, we derive the ground state energy of the highfrequency component and combine the results into an effective single-mode low-energy Hamiltonian.

Three-mode Hamiltonian and basis transformation

We start with the circuit Hamiltonian [START_REF] Smith | Magnifying quantum phase fluctuations with Cooper-pair pairing[END_REF]

Ĥ = 2E C 1 -ε 2 ( N 2 Σ + N 2 ∆ -2ε NΣ N∆ ) + 4ϵ C M 2 + 1 2 E L θ2 + ϵ L θ -φΣ -φ ext -1 2 θ ext 2 + ϵ L φ∆ -1 2 θ ext 2 -2E J cos φΣ cos φ∆ + 2εE J sin φΣ sin , (S6)
we have introduced symmetric and antisymmetric variables φΣ = 1 2 ( φ1 + φ2 ) and φ∆ = 1 2 ( φ1 -φ2 ). Note that N = NΣ N∆ M T are the conjugate Cooper-pair numbers to φ = φΣ φ∆ θ T . Unfortunately, model reduction using the Born-Oppenheimer approximation would be inefficient in this basis due to the strong coupling between θ and φΣ arising from the large value of ϵ L relative to the other energy scales.

To remedy this, we move to a basis where the dominant terms in the Hamiltonian are diagonal. The Hamiltonian in Eq. (S6) may be decomposed into two parts: a linear Hamiltonian that accounts for the energy of the capacitances and internal KITE inductance (proportional to E C , ϵ C , and ϵ L ), and a nonlinear potential (proportional to E L and E J ). In our parameter regime, we can regard the latter as a perturbation and proceed by diagonalizing the former, which has the form lim

E J →0 E L →0 Ĥ = 4 N T E C N + 1 2 φT E L φ E C = ϵ C   z -εz 0 -εz z 0 0 0 1   E L = 2ϵ L   1 0 -1 0 1 0 -1 0 1  
upon a gauge transformation to shift φΣ → φΣ -φ ext -1 2 θ ext and φ∆ → φ∆ + 1 2 θ ext . We have also defined

z = E C 2ϵ C (1-ε 2 ) = E C 2ϵ C + O(ε 2
) to parameterize the degree of capacitive loading, which diverges in the limit of vanishing junction capacitance. This linear Hamiltonian can be diagonalized with the transformation

φ = Λ φ Λ =   1 (z -1)ε -z 0 1 z(1 + z)ε 1 -ε 1   + O(ε 2 ) ,
which has the effect of shuttling the coupling from terms of order ϵ L to terms of smaller order E J and E L . We denote the transformed superconducting phases φ = φ θ ζ T . Note that in the limit of vanishing junction asymmetry and capacitance, z diverges and θ → φ. Recalling that the charging and inductive energy matrices transform according to

E C → Λ -1 E C (Λ T ) -1 and E L → Λ T E L Λ, we arrive at the Hamiltonian Ĥ = 4ϵ C 1 + z -1 N 2 φ + 2E C N 2 θ + 4ϵ C 1 + z N 2 ζ + ϵ L θ2 + ϵ L (1 + z) 2 ζ2 + 1 2 E L φ + ζ -ε θ + φ ext + 1 2 θ ext 2 -2E J cos φ -z ζ + (z -1)ε θ × cos θ + z(1 + z)ε ζ + 1 2 θ ext + 2εE J sin φ -z ζ sin θ + 1 2 θ ext + ϵ L O(ε 2 ) , ( S7
)
where we have introduced the transformed conjugate Cooper-pair numbers Nφ , Nθ , and Nζ ; and shifted φ → φ + φ ext + 1 2 θ ext . Note that the charging energy of the φ mode approaches ϵ C when the junction capacitance vanishes. On the other hand, when the shunting capacitance and junction asymmetry vanish, the Josephson potential terms in Eq. (S7) are identical to those in Eq. (S6).

Born-Oppenheimer approximation

In our parameter regime, ϵ L is the dominant energy scale. Therefore, the Hamiltonian in Eq. (S7) describes two high-frequency modes θ and ζ, and one low-frequency mode φ, which are all perturbatively coupled through the terms proportional to E J and E L . Our goal is to derive an effective Hamiltonian for the low-frequency φ mode. A method that is well suited for this purpose is the Born-Oppenheimer approximation, where one sets φ to a classical parameter φ in Eq. (S7) and removes the kinetic term proportional to N 2 φ , yielding the Hamiltonian Ĥθζ . One then needs to calculate the ground state energy E 0 (φ) of Ĥθζ . The Born-Oppenheimer approximation states that the effective potential of the φ mode is well captured by E 0 ( φ). However, an analytic expression for E 0 (φ) is unknown. Instead, we seek a low-order expansion of Eq. (S7) in order to analytically compute an approximate expression for E 0 (φ).

We recall that we are placed in the parameter regime where E L ≪ E J ≲ ϵ L . In the following, we will assume E J /ϵ L = O(λ) and E L /E J = O(λ), where λ is a small parameter. In the semiclassical limit of vanishing quantum fluctuations of θ and ζ, the effective potential of φ coincides with the full potential Û from Eq. (S7) evaluated at the equilibrium values θ0 and ζ0 . The equilibrium values themselves are found by solving ∂ Û ∂ θ θ0, ζ0 = 0 and

∂ Û ∂ ζ θ0, ζ0
= 0. Retaining only the lowest-order terms in λ and ε, we find θ0 = -

E J ϵ L sin 1 2 θ ext cos φ + O(λ, ε) ζ0 = E J z ϵ L (1 + z) 2 cos 1 2 θ ext sin φ + O(λ, ε) . (S8)
The low-frequency Hamiltonian can then be calculated-in the semiclassical limit of vanishing fluctuations of the phases θ and ζ-to be

Ĥφ = 4ϵ C 1 + z -1 N 2 φ + Û θ0, ζ0 = 2E C 1 + z N 2 φ + 1 2 E L φ + φ ext + 1 2 θ ext 2 - E 2 J 2ϵ L sin 2 1 2 θ ext - z 2 (1 + z) 2 cos 2 1 2 θ ext cos 2 φ -2E J cos 1 2 θ ext cos φ + 2εE J sin 1 2 θ ext sin φ + λE J O(λ, ε) + E C O(ε 2 ) . (S9)
Finally, shifting φ → φ-φ ext -1 2 θ ext and setting θ ext = π yields the single-mode model referenced in the main text,

Ĥφ π = 4ϵ C 1 + z -1 N 2 φ + 1 2 E L φ2 + E 2 J 2ϵ L cos [2( φ -φ ext )] -2εE J cos( φ -φ ext ) + λE J O(λ, ε) + E C O(ε 2 ) (S10)
after neglecting constant terms. We see that the shunting capacitance renormalizes the charging energy, while the inductive energy is completely unchanged. The first, second, and third terms in the potential represent the contributions of the inductive shunt, two-Cooper-pair tunneling, and asymmetric dc-SQUID, respectively. The Hamiltonian in Eq. (S10) is the main result of this section, and it explains both the approximate values of the effective parameters Ω, E J , δE J , and η of the low-energy Hamiltonian as well as the functional form. Overall, we sought a lowest-order result in the the parameters ε and λ, and completely neglected fluctuations of the high frequency modes. The parameters achieved in the experiment are not especially deep in this parameter regime-for instance E J /ϵ L = 0.58 and the fluctuations of θ are of order E C 2ϵ L 1/4 = 0.82 (similarly for ζ)-making the effective parameters in Eq. (S10) differ slightly from those found using a fit (see Tab. S2). Improved accuracy is likely possible using higher order corrections in λ and including fluctuations. 3b) at probe frequency 5.055 32 GHz. The phase response is periodic in two directions, which we identify with white arrows to be the independent external fluxes φext and θext.

EXTERNAL FLUX DEPENDENCE AND COHERENCE TIMES

Due to their layout, our physical flux bias lines (see Fig. 3b) inevitably couple to both circuit external fluxes φ ext and θ ext . For initial characterization, we probe the readout resonator at the single frequency 5.055 32 GHz and measure the reflected phase as a function of both bias currents. As the resulting pattern is 2D-periodic with many flux quanta visible, we can easily find the affine transformation that maps bias currents to φ ext and θ ext , as shown in Fig. S9. Additionally, we observe external flux drifts on the order of 1-2 % of a flux quantum on daily timescales, which we calibrate regularly using the procedure described in Ref. [START_REF] Smith | Magnifying quantum phase fluctuations with Cooper-pair pairing[END_REF].

In addition to the spectroscopy data shown in the main text, we acquire coherence data of the lowest frequency transition-the qubit-at the four external flux sweet spots, as shown in Fig. S10. As visible in Fig. 3, the qubit frequency varies by an order of magnitude from 636 MHz to 6.882 GHz as the external flux point is stepped from (θ ext , φ ext ) = (0, π) to (π, π) to (π, 0) to (0, 0); in other words, counter-clockwise around the plaquette in Fig. S10 starting at the green point. As the qubit frequency increases, the observed relaxation times decrease from 21.3 µs to 4.8 µs, indicating a roughly constant quality factor loss channel such as dielectric loss. The coherence times measured with a Ramsey sequence range from 0.8-6.6 µs and improve by about a factor of two using a single echo pulse (not shown). The dephasing is likely due to a combination of second-order flux noise, phase-slips in the array junctions, and photon shot noise in the readout resonator.

FIG. 1 .

 1 FIG.1. Principle of high-order photon processes. (a) Electrical circuit depicting a superconducting LC oscillator (black) shunted by a generalized Josephson junction (grey) that only permits Cooper-pair tunneling in pairs. The circuit is threaded by an external flux denoted ϕext. (b) Potential and energy levels (solid lines) of this nonlinear oscillator [Eq.[START_REF] Chang | Quantum nonlinear optics -photon by photon[END_REF] with parameters in Tab. I] and its linear equivalent (dashed lines) as a function of the superconducting phase difference at ϕext = π. Crucially, the adjacent transition frequency shift δn from the bare frequency ω0 alternates in sign when ascending the ladder. (c) Photon process diagrams for the first four interaction orders and corresponding interaction energies Jn.

FIG. 3 .

 3 FIG. 3. Experimental implementation. (a) Electrical circuit for the KITE (green) shunted by a superinductance (blue) and coupled through a shared inductance (purple) to an LC oscillator used for dispersive readout. The two small junctions have slightly different Josephson energies E ± J = (1 ± ε)EJ and charging energies E ± C = EC /(1 ± ε), with ε ≪ 1, due to junction fabrication variation. (b) Optical micrograph of the physical device, with false color indicating the constituent Josephson junctions. Aluminum and niobium electrodes appear in white and grey, respectively. (c)-(f) Two-tone reflection spectroscopy (background subtracted) of the four lowest transitions from the ground state along the four edges of the primitive cell in the two-dimensional external flux landscape (inset diagrams), and (g)-(j) accompanying readout spectroscopy. Theoretical transition frequencies (dashed lines) are obtained from numerical diagonalization of the three-mode circuit model with six fitted Hamiltonian parameters. The dominance of two-Cooper-pair tunneling along the θext = π edge strikingly suppresses the residual flux dispersion of the first four transition frequencies.

4

 4 is half the Kerr shift per photon and J n≥3 can be neglected. Hence ω n ≈ ω 0 + 2nJ 2 /ℏ and the transition frequency monotonically shifts for each added photon [Figs.2(c) and 2(e)]. This is in stark contrast with the regime of extreme phase fluctuations explored in this work where the transition frequency shift may alternate in sign for each added photon [Figs.2(d) and 2(f)]

  ) and3(b). It consists of a high impedance LC oscillator. The inductance, which we aim to maximize, is formed by a chain of 90 Josephson junctions, resulting in an inductive energy E L /h = 0.53 GHz. The capacitance, which we aim to minimize, originates EJ /h EC /h EL/h ϵL/h ϵC /h ε 3.92 6.23 0.53 6.73 5.65 0.032 Ω/2π EJ /h δEJ /h η 2.86 0.79 0.27 3.40 TABLE I. Extracted model parameters. The first line of parameters are found by fitting the spectral lines inFig. 3(c)-(f) to the three-mode circuit Hamiltonian of Fig. 3(a) in the absence of the readout mode [26]. The second line of parameters are found by fitting the measured transition frequencies at θext = π [Fig. 4(a)] to the effective one-mode Hamiltonian in Eq. (4). All energy scales are given in gigahertz.

  (c)] by inverting Eq. (3) [31]. Notably, we find that |J 2 | ≈ |J 3 |, which is consistent with the extracted η = 3.4.

FIG. 4 .

 4 FIG. 4. Spectral interlacing. (a) Adjacent transition frequencies obtained from two-tone measurements (open circles) and numerical diagonalization of the one-mode Hamiltonian (solid curves) along θext = π. (b) Transition frequency shifts δn and (c) interaction strengths Jn extracted from measurements (open circles) and analytic expressions depending on the fitted circuit parameters (fitted circles) at the starred external flux value. The transition frequency shifts alternate in sign while remaining much smaller than the transition frequency itself, directly corresponding to similarly-sized Hamiltonian coefficients for four-wave, six-wave, and eight-wave mixing.

  FIG. S5. Optical micrograph of the device and scanning electron microscope images of the three Josephson junction types from a nominally identical sample. (Green) one of the two small KITE junctions. (Blue) 14 of the array junctions that form both the internal KITE inductance and the inductive shunt. (Purple) 12 of the array junctions that form the shared inductance between the circuit and the readout resonator, as well as the self-inductance of the readout. All junctions are fabricated in one step using Dolan bridges.

  ) are Fock states |n⟩ where n is an integer. Using the formula: for n ≥ 0, â † |n⟩ = √ n + 1 |n + 1⟩ and for n ≥ 1, â |n⟩ = √ n |n -1⟩ and â |0⟩ = 0, their associated eigenvalue E n take the form

FIG. S8 .

 S8 FIG. S8.Comparison of models of the device. (a-b) Transition frequencies from the ground state obtained from numerical diagonalization of the three-mode Hamiltonian in Eq. (S6) and the one-mode Hamiltonian in Eq. (S10) with parameters in the first row of Tab. I, in addition to analytic results from RWA with the same parameters. Note the difference between these parameters and those obtained from a Hamiltonian fit (first and second rows of Tab. S2, respectively). (c) Adjacent transition frequencies emphasizing the differences between the three models.

FIG. S9 .

 S9 FIG.S9. Reflected readout resonator phase as a function of dc current in both loop and KITE flux bias lines (respectively the top and left current sources in Fig.3b) at probe frequency 5.055 32 GHz. The phase response is periodic in two directions, which we identify with white arrows to be the independent external fluxes φext and θext.

Acknowledgments: W.C.S acknowledges fruitful discussions with Agustin Di Paolo. The devices were fabricated within the consortium Salle Blanche Paris Centre. We thank Jean-Loup Smirr and the Collège de France for providing nano-fabrication facilities. This work was supported by the QuantERA grant QuCOS, by ANR 19-QUAN-0006-04. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement no. 851740). This work has been funded by the French grants ANR-22-PETQ-0003 and ANR-22-PETQ-0006 under the "France 2030 Plan."

Supplementary Material

DEVICE FABRICATION