Multimode physics of the unimon circuit

.


I. INTRODUCTION
Superconducting circuits are one of the most promising platforms to enable fault-tolerant quantum computing [1].However, reaching the stage where such devices become useful in practical applications still seems a major challenge, calling for gate fidelities and coherence times beyond the current state-of-the-art qubits, such as transmons [2][3][4].Even to achieve useful quantum advantage in the on-going noisy-intermediate-scale-quantum (NISQ) era [5], gate fidelities exceeding 99.99% for both singlequbit and two-qubit gates may be required, which has not been achieved in superconducting circuits yet.Thus, increasing the quality of superconducting qubits through design and fabrication is one of the greatest on-going technical challenges in the field.
Recently, different types of unconventional qubits combining desired features have been proposed as alternatives to transmons .One of them, the fluxonium, has demonstrated coherence times of the order of milliseconds [10], and averaged gate fidelities exceeding 99.99% for single-qubit gates [10] and 99.7% for two-qubit gates [8,11], thus providing alternative routes towards large-scale quantum computers.However, the involved architecture of the fluxonium may limit its reproducibility in fabrication.In addition, its low frequency requires special techniques for control, readout, and reset [9].
The unimon is another unconventional qubit recently proposed and tested experimentally [27].The unimon circuit exhibits simplicity since it consists of an inductively shunted single Josephson junction that can be biased by external flux.In contrast to the fluxonium, the unimon operates in the regime where the inductive energy of the shunt is mostly cancelled by the quadratic component of the Josephson potential.Consequently, the unimon circuit promotes not only high anharmonicity but also full insensitivity to dc charge noise and reduced sensitivity to homogeneous flux noise.
Despite the promising first experimental results of the unimon exhibiting 20% anharmonicities exceeding 0.7 GHz and 99.9% single-qubit gate fidelity [27], there is an urgent demand for a detailed understanding of its underlying physics.In particular, a comprehensive grasp of the fundamental multimode effects originating from the coplanar waveguide (CPW) resonator is important, since theoretical models employed in the description of the unimon to date have been restricted to the first and second-lowest normal modes of the system [27].Even though they provide a good qualitative agreement with the experimental results, these models do not fully capture the influence of the high-frequency modes.Hence, there is a demand in development of more involved models of the unimon.
In this paper, we develop the theory of multimode unimon circuits and address the physical phenomena induced by its high-frequency modes.Starting from the continuous distributed-element circuit, we quantize the system obtaining an auxiliary-mode Hamiltonian that is equivalent to the effective Hamiltonian obtained using the path-integral-based approach in Ref. [27].We then proceed with a partial linearization procedure [28][29][30][31][32][33] and find a renormalization of the Josephson energy [33,34] that was overlooked in Ref. [27].The Hamiltonian is then represented in a single-mode unimon basis that allows for efficient numerical diagonalization within the low-energy subspace.We put forward an efficient method to obtain corrections of the qubit energies and anharmonicities induced by the coupling of the qubit mode to a several unoccupied high-frequency modes.For typical unimon parameters, we find anharmonicity reductions of roughly FIG. 1.(a) Schematic representation of a unimon circuit featuring a single Josephson junction (EJ, CJ) that is displaced from the center of the coplanar waveguide (CPW) resonator (xJ ̸ = 0) with a total length of 2l.The CPW capacitance and inductance per unit length are denoted as C l and L l .The white regions symbolize silicon (Si), and the light blue areas represent niobium (Nb), excluding the Josephson junction composed of aluminum.The external magnetic flux biases through the first and second loops are labeled as Φext,1 and Φext,2.Voltage profiles for the first three modes (m = 1, 2, 3) are illustrated with solid lines.The displacement of the Josephson junction renders its coupling to all modes to be more evenly distributed, as emphasized by the discontinuity in the voltage profile.Please note that the color-coding for the modes remains constant throughout the figure.(b) The effective charging energy EC,m = c 2 m e 2 /(2C eff ), (c) transition frequency (Em,1 − Em,0)/h of mode m, (d) anharmonicity [see Eq. (15)] of mode m, and (e) cross-Kerr interaction [see Eq. (16)] between modes m − n as functions of the displacement xJ of the Josephson junction.For panels (b)-(e), the simulations were conducted using the parameter values from Table I with the number of low-lying modes set at M0 = 8.The results depicted in panels (c)-(e) are derived from the energy levels featured in Fig. 3(a).As a consequence, the avoided crossings lead to visible discontinuities in these results.
30% when up to eight modes are taken into account.
Previously, unimon has been studied in the special case where junction is located at the center of the circuit, leading to half of the modes being decoupled from the junction.Importantly, our results indicate that an asymmetric arrangement where the junction is offset from the center provides a rich profile of distributed nonlinearity, which is expressed by significant self-and cross-Kerr couplings between the modes.Particularly, the numerical findings from both two-and three-mode scenarios display moderate nonlinearity across all modes concurrently, hinting at the potential for multimon-like qubit operations [35][36][37] within unimon circuits.However, to achieve high-fidelity multiqubit operations, unimon circuits with more intricate designs or additional components are likely needed.We further study the accuracy of our numerical findings by analytically solving for self-and cross-Kerr interactions within the harmonic-oscillator basis.
This paper is organised as follows.In Sec.II, we introduce our model and notation, detailing the linearization procedure and obtaining spectrum of the qubit.In Sec.III, we describe the effects of the high-frequency modes on the unimon-qubit mode comparing them with the results predicted by the single-mode approximation and with the auxiliary-mode models introduced in Ref. [27].In Sec.IV, we study the multimode structure of the unimon, focusing on the self and cross-Kerr terms.Our conclusions are presented in Sec.V.

II. MULTIMODE MODEL FOR THE UNIMON
We study a single-junction unimon circuit with multiple modes taken into account.The system is schematically illustrated in Fig. 1(a).Our primary motivation lies in discerning how the position of the Josephson junction within the half-wavelength resonator influences the nonlinearity of the modes.Key indicators of this nonlinearity, namely the mode effective charging energy, transition frequency, anharmonicity, and cross-Kerr coupling between modes, are presented in Fig. 1(b)-(e).To accurately compute these values, we introduce an effective theoretical model, which is elaborated upon in the subsequent sections.

A. Effective multimode unimon Hamiltonian
As a starting point for our multimode treatment of the unimon circuit, we use a Hamiltonian that comprises a nonlinear mode with Ψ describing the magnetic flux difference across the Josephson junction and being conjugate to the charge Q, and M included linear auxiliary modes with fluxes χm that are conjugate to charges Ξm , where m = 1, ..., M .The relevant nonzero single-operator commutation relations these operators satisfy are [ Ψ, Q] = iℏ and [ χk , Ξm ] = iℏδ km .A detailed derivation for the Hamiltonian is provided in the Appendix A, and hence, we begin our treatment from the auxiliary-mode Hamiltonian where C eff is the effective capacitance, L l is the resonator inductance per unit length, 2l is the resonator length, and the auxiliary-mode coupling strengths assume the form where v = 1/ √ L l C l is the phase velocity, C J is the capacitance of the Josephson junction, C l denotes the capacitance per unit length, and x J ∈ [−l, l] represents the location of the Josephson junction in the CPW resonator.
The auxiliary-mode Hamiltonian in Eq. ( 1) contains M auxiliary modes.Although this expression becomes exact as M → ∞, numerical calculations necessitate limiting the number of modes at a finite M .How to choose M is a relevant question for our multimode model and is addressed in more detail in Sec.III.In addition, we note that the version of the auxiliary-mode Hamiltonian introduced here differs slightly from the one presented in Ref. [27].The distinction lies in the dc flux offset parameter ϕ 0 that we introduce (Appendix A) to facilitate our derivation of the multimode model.
The auxiliary modes couple to the nonlinear mode inductively, providing corrections to the unimon qubit energy levels.The auxiliary-mode Hamiltonian in Eq. ( 1) is well-adapted for studying the unimon circuit if our interest is focused only on the lowest mode, as in Ref. [27].Since the auxiliary-mode frequencies are integer multiples of the lowest auxiliary mode, Ω m = πmv/(2l), only the few lowest modes are energetically close enough to the nonlinear mode Ψ to significantly interact with it.Furthermore, the location of the Josephson junction can be chosen strategically to leave certain auxiliary modes uncoupled, thus easing the load for numerical calculations.In an optimal scenario, this approach allows for fairly accurate results for unimon qubit energy levels to be obtained by incorporating just a single auxiliary mode, effectively reducing the problem to solving a twodimensional Schrödinger equation in the flux basis [27].However, extending the consideration from the energy levels of the lowest mode to those of a multimode system necessitates an alternative approach, primarily due to the rapid escalation of computational demand as more auxiliary modes are included.
To adapt the Hamiltonian of the unimon circuit for numerical analysis in cases involving multiple modes, we divide the auxiliary-mode Hamiltonian in Eq. ( 1) into linear and nonlinear parts, Ĥaux = Ĥlin + Ĥnl , by expanding the nonlinear Josephson term and moving the resulting quadratic term to the linear part (Appendix B).To elucidate different modes in the system and simplify subsequent analysis, we find the classical normal modes of the linear part, Ĥlin , using a basis transformation.This process, which effectively removes the linear coupling between the flux operators in Eq. ( 1), is detailed in the Appendix B. We note that similar linearization procedures to find the normal modes of system have been employed in earlier works, as seen in Refs.[28][29][30][31][32][33].
The basis transformation gives rise to the normal-mode flux operators and the corresponding conjugate charge operators which are denoted by φm and qm , respectively.These new operators continue to satisfy the canonical commutation relations [ φk , qm ] = iℏδ km .The transformation also yields the normal-mode representation of the magnetic flux across the Josephson junction: where the constant factors c m are coefficients determined by the transformation and describe the contribution from each normal mode to the overall magnetic flux.In addition, the diagonalization process uncovers the normalmode frequencies, denoted by ω m /(2π), m = 1, ..., M + 1. Insertion of Eq. ( 4) into the nonlinear part Ĥnl gives us the full normal-mode representation of the auxiliary-mode Hamiltonian as where we have defined the Josephson inductance to represent the effective inductance of the m:th normal mode.Above, we described the Hamiltonian of the unimon circuit in the normal-mode basis of the linearized version of the circuit.However, a challenge still remains, since finding the energy levels of the unimon qubit, while including the interactions from the higher number of modes, requires solving a high-dimensional Schrödinger equation, and thus, an improved basis to find the solution is needed.
In order to simplify the notation, we introduce dimensionless operators nm = (q m /c m )/(2e) and φm = 2πc m φm /Φ 0 , where e is the elementary charge.We divide the Hamiltonian in Eq. ( 5) into parts describing the singlemode unimon Hamiltonians Ĥm and the interaction part Ĥint , leading to These constituent Hamiltonians can be expressed as and where we have defined an effective charging energy of the m:th mode as E C,m (φ 0 ) = c 2 m (φ 0 )e 2 /(2C eff ), effective inductive energy as E L,m (φ 0 ) = Φ 2 0 /(2π) 2 /[ Lm c 2 m (φ 0 )], and the dc phase across the junction is denoted as φ 0 = 2πϕ 0 /Φ 0 .Note that the coefficients c m introduce the dependence of E C,m and E L,m on φ 0 .This normal-mode representation of the auxiliary-mode Hamiltonian is beneficial for several reasons.On one hand, it fully separates the single-mode components from the interaction terms.This becomes particularly clear on the first row of Eq. ( 9), where the first term cancels the concealed single-mode terms in the second term (see Appendix C).On the other hand, the single-mode unimon Hamiltonian Ĥm can be diagonalized efficiently using the one-dimensional Schrödinger equation, Ĥm |j m ⟩ = E m,j |j m ⟩, where |j m ⟩ and E m,j denote the j:th eigenstate of the m:th mode and the corresponding eigenenergy.Moreover, the eigenstates of Ĥm form a set of basis states that can be used to represent the interaction part Ĥint in a matrix form, that subsequently allows an efficient diagonalization of the Hamiltonian in Eq. ( 7) for the purposes of studying the effect of the other modes on the energy levels of the unimon qubit.In practice, this is accomplished by applying the generalized trigonometric sum relation cos for a cosine of sums, and subsequently calculating matrix elements of three different types, j m cos φm j ′ m , j m sin φm j ′ m , and ⟨j m | φm |j ′ m ⟩, in Eq. ( 9).To avoid any possible confusion, we emphasize that the last summation in Eq. ( 10) represents a sum over all possible subsets of {1, ..., N } with a size of k.The product over i / ∈ A refers to all elements in {1, ..., N } not included in A. When k = 0, A is an empty set, and the product over i / ∈ A encompasses all elements of {1, ..., N }.

B. Energy cutoff for Hilbert space
After establishing a matrix representation in the singlemode unimon basis, we need to manage the dimensions of the Hilbert space before proceeding with an efficient diagonalization of the matrix.To this end, we restrict the total size of the Hilbert space, which ensures computational feasibility and accuracy of our numerical approximations.
Although the number of energy levels for each mode is theoretically infinite, ideal quantum computation takes place in a finite-dimensional space and, therefore, is compatible with the concept of an energy cutoff E cutoff .Consequently, we disregard all eigenstates of the single-mode Hamiltonian Ĥm with energies exceeding E cutoff .The value of the energy cutoff is determined by the convergence of the low-energy eigenstates of the full Hamiltonian that we aim to accurately model.Since the bare frequencies of the modes increases with the mode number m, we only need to consider states beyond the ground state for a limited number of modes.This is attributed to the fact that the energy required to excite such a mode exceeds the established energy cutoff.
We categorize the complete set of modes (m = 1, 2, ..., M + 1) into two distinct groups, referred to as the lower modes and higher modes.The lower modes, defined by an integer M 0 and m ∈ {1, ..., M 0 }, include all modes where any excited states are considered.Conversely, for the higher modes, for which m ∈ M 0 + 1, ..., M + 1, only their vacuum state is included due to the energy cutoff.
To simplify numerical computations, we assume that the system is operated at a flux sweet spot where Φ diff = Φ 0 /2 and φ 0 = π.This condition results in a symmetry ⟨φ and hence the expectation values of operators anti-symmetric in φm vanish for the vacuum sate.For example, 0 m sin φm 0 m = 0 and ⟨0 m | φm |0 m ⟩ = 0. Interestingly, this implies that the impact of the higher modes on the system primarily contributes to a renormalization of the Josephson energy, denoted as A more comprehensive exploration of this renormalization effect and its implications can be found in Sec.III.
After incorporating all of the above-described steps, we arrive at the final form of the total Hamiltonian where we have used Eq. ( 10) and truncated the Hilbert space based on the energy cutoff.We consider this expression to be one of the main results of this paper.
In our model, the Hamiltonian described by Eq. ( 12) is expressed in a matrix form using the single-mode unimon basis, where matrix elements are expressed as and each state |i 0 j 1 ...k M0 ⟩ must meet the energy cutoff condition By carefully selecting the energy cutoff, the diagonalization of the matrix can be accomplished with adequate numerical efficiency and accuracy.The above-introduced process of solving the eigenstates of the unimon circuit is visually summarized in Fig. 2.

C. Labeling of eigenstates
Upon diagonalizing the matrix, we acquire a new set of eigenstates along with corresponding energy eigenvalues, which can be interpreted as perturbed versions of the single-mode unimon eigenstates.To enhance our understanding of the effects of the interactions between the modes, it is fruitful to study the energy differences between the perturbed and non-interacting scenarios.
Although visual inspection of the energy levels can yield insight in specific cases, this approach tends to become increasingly demanding in general.To streamline  maximum overlap with |α n ⟩, and label |α n ⟩ as the perturbed counterpart of |i 0 j 1 ...k M0 ⟩.Note that our labeling method may produce ambiguous results.This ambiguity arises from transverse-type interactions introduced by the interaction term, leading to strong hybridization between eigenstates near suitable degeneracy points.Figure 3(a) illustrates the effects of the hybridization in the energy level diagram, exhibiting avoided crossings between levels.Where eigenstate labeling is applied on quantities such as anharmonicity or cross-Kerr interaction, the effects of hybridization are exposed through sudden discontinuities as found in Fig. 1(d)-(e).

III. MULTIMODE EFFECTS IN THE UNIMON CIRCUIT
Our next step is to analyze how the modes beyond the lowest mode affect the energy levels and anharmonicity of the unimon qubit.First, we focus on the anharmonicity calculated by using the renormalization model.In contrast to the previous single-mode approximations, as shown in Eq. ( 8), here we use the renormalized Josephson energy, denoted as ẼJ .In addition, we discuss the anharmonicity results from the multimode model detailed in Sec.II and juxtapose them with the outcomes of the renormalization model.As a further point of comparison, we also discuss the distinctions between these findings and those acquired using the auxiliary-mode model of Ref. [27].

A. Effect of renormalization
As indicated in Eq. ( 11), at the sweet spot φ 0 = π, the higher modes influence the energy levels of the unimon even in their vacuum state, where these modes still exhibit zero-point phase fluctuations.These fluctuations couple to the unimon mode, effectively causing a renormalization of the Josephson energy within the system [33,34].Given that the renormalization coefficient of Eq. ( 11) is simply a product of the vacuum state expectation values of the cosine function, it follows that ẼJ /E J ≤ 1.As a result of the renormalization, we thus observe a decrease in the anharmonicity of the unimon qubit.
Interestingly, this decrease in anharmonicity due to the renormalization of the Josephson energy offers an explanation for the discrepancies observed between the two models used in Ref. [27]: the single-mode approximation and the auxiliary-mode model based on path integrals.The auxiliary-mode model in Eq. ( 1), which includes the nonlinear mode and a few lowest linear modes, is expected to provide a more accurate results than the single-mode approximation model, which only incorporates one normal mode [Eq.( 8)].
In Fig. 4(a), we display the discrepancy in qubit anharmonicities.Between the single-mode and auxiliary-mode (M = 2) approximations, the difference is roughly 20% at x J = 0.By including the renormalization, we observe an approximate 20% decrease in the anharmonicity compared to the single-mode approximation, bringing the result appreciably close to the predictions of the auxiliary-mode approach.In Fig. 4(b), we further detail the behavior of the renormalization coefficient, which changes only weakly as a function of x J .
It is worth mentioning that based on physically motivated mode-cutoff frequncies [33,38], the number of modes M included in the renormalization process should be finite.Here we introduce a cutoff that is based on the magnitude of the superconductor gap parameter ∆ gap .(11)].These numerical calculations are carried out using the parameters listed in Table I.
Namely, we set Ω m ≤ 2∆ gap , which imposes a limit on the number of modes [33,39,40].However, the renormalization coefficient appears to be insensitive to the precise number of modes included and seems to converge as M → ∞.This is evidenced by the inset of Fig. 4(b), where the contribution from each mode around the cutoff frequency (M ≈ 100) is negligible.More detailed discussion regarding the convergence of the renormalization coefficient is found in Appendix D.
Our numerical findings shown in Fig. 4 demonstrate that the renormalization of the Josephson energy leads to anharmonicities close to those obtained from the flux basis solution of the auxiliary-mode Hamiltonian.Importantly, compared to the single-mode-approximation model, the renormalization model achieves a comparable level of accuracy to the auxiliary-mode approach (M = 2), with virtually no increase in computational load.

B. Effects beyond renormalization
In addition to the renormalization effect of vacuum states, we aim to understand how the energy levels of the unimon qubit are influenced when excited eigenstates in modes m > 1 are included within the computational Hilbert space after truncation.Therefore, we employ the comprehensive multimode method to solve for the energy levels and anharmonicities of the full Hamiltonian in the single-mode basis as shown in Eq. (13).Our findings show that the anharmonicity decreases further in comparison to the renormalization model.This indicates that applying the multimode model is important for obtaining accurate results.
We also propose that this model exhibits higher numerical accuracy with a given number of computational resources compared to solving the auxiliary-mode model in the flux basis.This claim is supported by two factors.First, solving the auxiliary-mode Hamiltonian becomes in general very challenging for increasing M > 2 owing to the exponential increase in the dimension of the computational Hilbert space.Second, as depicted in Fig. 4, the results derived from the auxiliary-mode model appear to converge towards those obtained from our multimode model with increasing M .However, due to the hefty computational demands imposed by the auxiliary-mode model, we are unable to verify this comprehensively.Note that, in special cases such as x J /l = 0, we are capable of solving the auxiliary-mode model for M = 4 since only auxiliary modes with an even ordinal number exhibit nonzero coupling with the nonlinear mode.In addition, we emphasize that the multimode model requires significantly less computational time than the auxiliary-mode model to achieve a comparable accuracy.More details on the accuracy of the multimode model can be found in Appendix E.

IV. INTERACTIONS BETWEEN THE LOWEST MODES OF THE UNIMON CIRCUIT A. Energy levels and avoided crossings
Let us leverage the model developed in Sec.II to its full extent by examining the interactions among the three lowest modes of the unimon circuit.This involves solving the multimode Hamiltonian in the unimon basis as given in Eqs.(13), and subsequently employing our labeling method to map the obtained eigenstates onto the states of the unimon basis.The energies of the labeled eigenstates are denoted as E ijk , where i, j, and k represent the number of excitations in the first, second, and third modes, respectively.In Fig. 3, we show all obtained energy levels as functions of x J , E J , and Z c .Here, Z c = L l /C l represents the characteristic impedance of the transmission line forming the resonator.Within the energy levels, we highlight the labeled states where each of the initial three modes has at most one excitation.In Fig. 3(a), the energy levels manifest fairly intricate interactions with each other as a function of x J .We observe a correlation between the first excited states and the effective charging energies shown in Fig. 1(b), which is evident in the alignment of the local extrema.The effective charging energy for the mode depends on the weight coefficient c m which is related to the coupling between the Josephson junction and mode m [see Eq. ( 4)].Figs.3(b)-(c) illustrate that an increase in either the Josephson energy or the characteristic impedance results in a decrease in energy for the first excited states.This behavior can be attributed to the cancellation of quadratic terms as either E J or Z c increases.As these quadratic terms decrease, the potential becomes less steep, which in turn narrows the energy gap between the ground state and the excited states.However, the energies of the excited states of the third mode remain largely unchanged, given their weak coupling with the Josephson junction.As discussed in Sec.II C, these interactions often manifest as avoided crossings.An interaction exemplifying this can be observed between energy levels E 010 and E 300 at x J /l = 0.42, visualized in the inset of Fig. 3(a).
In qubit operations, avoided crossings are undesirable.If qubit energy levels become hybridized, there is a risk of unintended leakage into other participating states.This can adversely affect the coherence time of the qubit.Consequently, it is preferable to operate in parameter regimes where hybridization of the qubit states is minimized.The prevalence of avoided crossings increases with energy, rendering the high-frequency normal modes (m ≥ 2) challenging for qubit operation.Nonetheless, in most parameter regimes, the qubit of the lowest mode (m = 1) remains unaffected by these avoided crossings up to its second excited state.The notation αm(i, j) corresponds to the anharmonicity of the m:th mode (m = 1, 2, 3), with the other two modes having i and j excitations.The left argument (i) pertains to the lower mode of the two arguments.The cross-Kerr interaction is denoted by Knm(i), where n and m (n, m = 1, 2, 3) represent the modes between which the interaction is computed.The argument (i) equals to the number of excitations in the mode not primarily involved in the interaction.The simulation utilize the parameter values of Table I.

B. Anharmonicities and cross-Kerr interactions
We employ the labeled energy levels as a basis to compute key properties, such as the anharmonicity and cross-Kerr interaction.For the lowest mode m = 1, the anharmonicity is defined as and the cross-Kerr interaction between the modes 1 and 2 is defined as In general, the subscripts in α m and K mn define which modes are considered to change their excitation number in analogy with the above equations.Although it is customary to consider the anharmonicity and cross-Kerr interaction to involve only one and two modes, respectively, higher-order terms become significant when excitations of  15) and ( 16), respectively, and in Fig. 5.For the the parameter values that are not swept, we use Table I, except for the Josephson junction location set at xJ/l = 0.62, as marked by the vertical dashed line in panels (a) and (d).
the other modes are allowed.In this context, the anharmonicity of a mode depends on the occupation number of the other two modes, while the cross-Kerr interaction between two given modes is influenced by the state of the third mode.We consider these effects in more detail in Sec.IV C.
By employing Eqs. ( 15) and ( 16), we compute the anharmonicity and cross-Kerr interaction for scenarios where only the two lowest modes are allowed a single excitation.Consequently, each mode exhibits two state-dependent anharmonicities and the cross-Kerr interaction between these two modes can be characterized by a single value.The corresponding results for the two-mode scenario are displayed in Fig. 5. Our observations indicate that within the two-mode framework, it is feasible to identify a set of parameters that distribute the nonlinearity relatively evenly across both modes, for example, those given in Table I.Under this parameter configuration, if the other mode is in the vacuum state, the anharmonicity for each mode remains above a 50-MHz threshold.Furthermore, the cross-Kerr interaction energy lies around 200 MHz, peaking approximately at the same value of x J as α 2 (0, 0).However, it is evident that the anharmonicity of the lowest mode consistently surpasses that of the second mode.Transferring an excitation to a mode appears to have a diminishing effect on the anharmonicity of the other mode.The relative difference, between the anharmonicities with and without excitations present in the other mode, continues to increase with either increasing Josephson energy or characteristic impedance, even though the overarching trend is an increase in anharmonicity.Intriguingly, the value of α 1 0) even begins to decrease once E J or Z c exceed a certain threshold.More details on this effect is revealed in Sec.IV C.
Our numerical calculations also extend to anharmonicities and cross-Kerr interactions for the three-mode scenario, the results of which are illustrated in Fig. 6.In our efforts to spread the nonlinearity across all three modes, we adjusted the location of the Josephson junction to x J = during the sweeps with respect to E J and Z c .
the case of three modes, the number of anharmonicity parameters per mode increases to four, owing different combinations of excitations in the other two modes.Similarly, the number of cross-Kerr interaction parameters for any pair of modes doubles, reflecting the influence exerted by the state of the remaining mode.The results highlight the challenge in discovering parameters that ensures all anharmonicities exceed a 10-MHz threshold.As discussed in the context of the two-mode case, augmenting E J or Z c is not a feasible solution as it would result in a decrease in the lowest anharmonicity and cross-Kerr interaction.This outcome is visible for the anharmonicity α 2 (1, 0) as shown in Fig. 6(b)-(c), and for the cross-Kerr interaction K 23 (1) as depicted in Fig. 6(e)-(f).

C. Analysis of Kerr-type terms
Although our model lends itself well to numerical simulations, it does not facilitate analytical derivations.This limitation stems from the fact that the eigenstates of the single-mode unimon Hamiltonian seem not analytically solvable in general.Therefore, to analytically investigate the anharmonicities and cross-Kerr interactions, we make use of the harmonic-oscillator basis [30]. Here, m=1 e −λ 2 m /2 represents the renormalized Josephson energy in the harmonic-oscillator basis [33,34], where λ m = 2 E C,m /(ℏω m ) signifies the zero-point fluctuations of the m:th mode.We refer to the different modes with the indices n, m, and k.Under these definitions, the analytical expressions for the anharmonicity and cross-Kerr interaction are derived as where N m denotes the occupation number of the m:th mode and the self-Kerr and cross-Kerr interaction param-eters are given by For more details on the derivation of these expressions, refer to Appendix G.

Modes in vacuum state
We begin our analysis with Eqs. ( 19) and ( 20) that represent the case without additional excitations in the system.To gain insight into the behavior of the self-Kerr interaction with respect to x J , we represent K mm in a more suggestive form, K mm ∝ (E C,m /ω m ) 2 , where both E C,m and ω m are dependent on x J as depicted in Fig. 1(b)-(c).
Out of these two quantities, the effective charging energy E C,m is notably more sensitive to variations in x J .As a result, E C,m largely determines the dependence of the anharmonicity of mode m on x J .Interestingly, the dependence of E C,m on x J follows the magnitude of the voltage discontinuity of mode m across the junction, suggesting that E C,m serves as an indicator of the coupling strength between mode m and the Josephson junction.
On the other hand, the transition frequency, ω m /(2π), significantly influences the level of the anharmonicity for each mode.A higher transition frequency corresponds to reduced anharmonicity, as depicted in Fig. 1(d).Note that we do not consider the dependence of E * J on x J in detail because it is relatively weak, as illustrated in Fig. 4(b).
The analytical form of the cross-Kerr interaction K nm reveals a simple, yet important relationship with the self-Kerr interactions K mm and K nn .Specifically, the cross-Kerr interaction K mn is twice the geometric mean of the self-Kerr interactions K mm and K nn .This indicates that the most efficient way to increase the cross-Kerr interaction between the modes, is to augment the self-Kerr interaction of both modes by an equal amount.Such conclusions regarding the behavior of the cross-Kerr interaction also applies for the parameter sweeps respect to E J and Z c .
We continue by examining the parameter sweep with respect to the Josephson energy E J .In this case, the self-Kerr term of mode m can be expressed as K mm ∝ ẼJ (E C,m /ω m ) 2 , highlighting the dependency of the self-Kerr on E J .From the Fig. 3(b), we observe that the transition energies from the vacuum state E 000 to the first excited states E 100 , E 010 , and E 001 decrease with increasing E J .Consequently, this leads to an increasing effect on K mm , a trend also evident in the numerical results.
In addition, the renormalized Josephson energy has an enhancing effect on K mm , even if the renormalization coefficient may decrease with increasing E J .As a result, the self-Kerr interaction strength is anticipated to increase more rapidly with respect to E J for the lowest than for the other modes, a pattern consistent with our numerical calculations.Another approach to reach this conclusion is by examining the quadratic part of the singlemode Hamiltonian in Eq. ( 8), given by (E L,m − E J ) φ2 m /2.As E J /E L,m → 1, this term approaches zero, while the magnitude of the nonlinear terms grows [27].
For the sweep with respect Z c , we expect similar behavior of the self-Kerr interaction to that for E J .This is substantiated by the observation that when the operators in the auxiliary-mode Hamiltonian of Eq. ( 1) are rescaled, the outcomes of increasing either E J or Z c are approximately equal in the linear approximation (Appendix F).At the sweet spot, the linear effects are pronounced because the main cause of the increase in self-Kerr interaction is the cancellation of quadratic flux term as E J /E L,m → 1.Nevertheless, higher-order terms indicate different behavior between the characteristic impedance Z c and the Josephson energy E J .Such disparities are particularly evident in the case of high-frequency modes, where the effect of the cancellation diminishes.For additional details, refer to Appendix F.

Effect of excitations
We turn our attention to the impact of excitations on both anharmonicity and cross-Kerr interactions.The origins of these effects can be traced back to nonlinear Hamiltonian terms, which are elaborated upon in Eqs.(17) and (18).
Introducing a single excitation into the k:th mode yields two notable outcomes.First, there is a negative correction to the anharmonicity of mode m.Second, a corresponding negative correction appears in the cross-Kerr interactions between modes m and n.The relative decline for both of these effects is captured mathematically as which reveals that the decrease caused by the excitation is at its maximum where K kk /E * J attains its peak value.Furthermore, the equation suggests that a mode with a larger self-Kerr interaction induces a larger relative correction.These effects are evident in Fig. 5 for the anharmonicity and in Fig. 6(d)-(f) for the cross-Kerr interaction.
Finally, we examine the anharmonicity when another excitation is introduced to the previously unoccupied mode n.The difference between this and the preceding case can be expressed as The magnitude of the zero-point fluctuations in mode k, which is the mode excited for Eq. ( 21), determines the sign of the correction.Although the strength of the zero-point fluctuations in mode n does not influence the sign, it acts as a scaling coefficient, influencing the magnitude of the correction.In light of this analysis, the observed trends in Fig. 6 align well with our expectations.For a majority of the parameter configurations, the lowest mode shows that α 1 (1, 1) is the smallest anharmonicity.For the other two modes, we consistently find that α m (1, 1) > α m (1, 0).From a qualitative perspective this makes sense since the lowest mode has the largest zero-point fluctuation.However, quantitatively this behavior is not explained by Eq. ( 22), since positive correction requires λ 2 1 > 1 which is not satisfied in this case.
It is important, however, to recognize the limits of this analytical methodology.Although instrumental in elucidating the qualitative behavior of Kerr-type interactions, the approach does not encompass transverse-type or other rotating-wave interactions.These interactions may be significant, especially in systems with notable nonlinearity.

V. CONCLUSIONS
In this work, we explored the effects arising from the multimode nature of the single-junction unimon circuit.Whereas our primary focus was on the impact of the high-lying modes on the lowest mode, we also investigated the influence of the nonlinearity on other modes for different locations of the Josephson junction.To facilitate our study, we developed a theory of multimode unimon circuits.Leveraging this framework, we determined the energy spectrum for several low-lying modes using numerical diagonalization in the low-energy subspace.
Our model markedly differs from the one presented in Ref. [41] that describes the multimode physics in Josephson-junction array fluxonium circuits.While the effectiveness of their model relies on the decoupling of the qubit mode from all other modes, the distinctiveness of our model arises from the large energy separation between the qubit subspace and the high-frequency modes, which is a consequence of the CPW structure.Furthermore, the methods of numerical diagonalization for fluxonium qubits presented in Refs.[14,32] differ from our approach, particularly with respect to the chosen basis.In Ref. [14], diagonalization is executed in the normal-mode flux basis, which enforces a mode cutoff that includes only two modes.On the other hand, the Hamiltonian in Ref. [32] is expressed in the harmonic-oscillator basis, leading to a less efficient convergence of the low-energy eigenstates when compared to our single-mode unimon basis.By adjusting the choice of basis to the problem at hand, we expect that the presented model is also applicable to a class of multimode qubits known as noise-protected qubits [22][23][24][42][43][44].
Our findings reveal that multimode effects introduce significant corrections to both the transition frequency and anharmonicity of the unimon qubit.Utilizing our comprehensive model, the decrease in anharmonicity can be as much as 30% when the junction is centrally lowith more pronounced effects observed at other positions.We also showed that by solely considering the vacuum states, the multimode effects can be condensed into a single coefficient, leading to a renormalization of the Josephson energy.The decrease of the anharmonicity owing to the renormalization can be compensated by the choice of an increased unnormalized Josephson energy.Compared to the single-mode approach, this method seamlessly incorporates certain multimode effects without greatly increasing the required computational resources.
We found that the second and third modes present a diverse nonlinearity profile, but their peaks do not coincide.This misalignment poses a challenge in identifying an operational regime where all three of the lowest modes exhibit substantial nonlinearity.Yet, for the first and second modes, we identified a position with notable nonlinearity in both.
Our work significantly expands upon the prior theoretical descriptions of the unimon [27] by introducing a theoretical framework to systematically investigate the multimode physics of the unimon circuit and its potential applications.An interesting application for unimon circuits is to encode multiple qubits into a single device in the spirit of multimon qubits [35][36][37].However, although we identified parameter configurations where two modes simultaneously exhibit significant nonlinearity, we did not find promising parameters for particularly highfidelity qubits.In our future research, we aim to apply our theoretical framework for more complex unimon circuits.For example, by implementing a ring-like geometry and adding more junctions, we expect to witness behavior similar to transmon-based multimons, characterized by a relatively uniform distribution of nonlinearity among the modes with high intrinsic anharmonicity.
We initiate our investigation by formulating the total Lagrangian for the continuous distributed-element circuit, presented as where we have defined which describes the flux of the grounded CPW resonator and Ψ(t) = ψ 2 (x J , t)−ψ 1 (x J , t) denotes the flux difference across the Josephson junction.The Lagrangian density L of the CPW and the Lagrangian L describing the flux across the junction are defined as where Φ diff is the total flux difference across the loops and the system obeys the boundary condition ψ(−l, t) = ψ(l, t) = 0.
To deduce the equations of motion, we employ the Euler-Lagrange equations, defined as where i ∈ {1, 2} and ψ J,i ≡ ψ i (x J , t).To utilize the Euler-Lagrange framework, we derive the subsequent equations of motion where the phase velocity is denoted as v = 1/ √ L l C l .

Variable elimination in frequency domain
Our interest is to eliminate the variable ψ i and express the equation of motion solely with the variable Ψ.In order to achieve this, we move to a frequency domain by using Fourier transformation.
d n K(ω)/(dω n ) ω=0 and the latter term in the second row of Eq. (A17) are zero.
In addition, since the only lumped element which depends time-derivatives of flux is capacitance, we may utilize the low-frequency assumption made with the pole expansion to approximate the convolution by neglecting time-derivatives ∂ n t Ψ(t) when n > 2. Interestingly, this low-frequency approximation made in the convolution expression in Eq. (A20) becomes exact in the limit of M → ∞.In this limit, Mittag-Leffler theorem [45] allows us to express the kernel in Eq. (A17) as where all contributions from kernel K(ω) to the timederivatives in the convolution term vanish.Similar technique has been previously utilized in the context of network synthesis of prescribed impedance functions [46].The approximation of convolution, as given in Eq. (A20), leads us to a more tractable expression for the full equation of motion in Eq. (A16) which takes the form where we used the definition for C eff given in the Eq. ( 2) of the main text and defined To handle the temporally nonlocal term, we employ a set of auxiliary modes, denoted as where ξ m = √ 2r m Ω m C eff .This helps us to express Eq.(A22) as For a complete description of the system dynamics, we derive the equations for each auxiliary mode.By utilizing Eq. (A24) and the convolution theorem, we obtain By introducing auxiliary modes, we have successfully eliminated the need to compute temporal convolutions.

Classical treatment of the dc flux
In the pursuit of computing the normal modes of the system, it is convenient to redefine the fluxes such that they vanish at the minima of their effective potentials.To this end, we define the shifted flux and the auxiliary variables as follows: Ψ(t) ′ = Ψ(t) − Φ diff + ϕ 0 represents the total flux as a sum of its dynamic component Ψ ′ (t) and a shift ϕ 0 .Similarly, for the auxiliary modes, we define χ ′ m (t) = χ m (t) + x m .Substituting these into our earlier equations, we obtain It is clear that the right side of these equations captures the time-independent behavior, with the left side containing time-dependent part of the system.
Conveniently, the time-independent parts of Eqs.(A27) and (A28) can be solved independently from the timedependent parts, giving us a set of new equations From these equations, by eliminating the variables x m , we obtain the relation This equation links the dc flux across our junction to the external magnetic flux, encapsulating the dc flux behavior of our system in the presence of an external magnetic influence.Note that the Eq.(A31) becomes multi-valued if 2lL l /L J > 1 [47][48][49], a parameter region we aim to avoid.
where the normal-mode frequencies are denoted as ω m /(2π).Furthermore, the normal-mode decompositions for operators Ψ and χm take the form and thus, the coefficients c m in Eq. ( 4) of the main text are given by c m = [U] 0,m−1 .
Appendix C: Division to single-mode and interaction parts Starting with the assumption that φ 0 = π and inspecting Eq. ( 10) of the main text, we note that only the term containing only cosines is a single-mode term.All other terms involve products of sines, which introduce coupling between the modes.We proceed by expressing the product of cosines as The significance of cos φm − 1 arises from the cancellation of the constant term.Consequently, the last summation term inevitably includes only interaction terms, whereas the second summation term holds the single-mode contribution.By applying the above steps to Eq. ( 9), the single-mode contributions conveniently cancel each other out, and hence only interaction terms remain.

Appendix D: Scaling of the renormalization coefficient
Here, we study the scaling properties of the renormalization coefficient ẼJ /E J in more detail.We put forward a convincing argument that the mode cutoff frequency used in the multimode model of the main text is reasonable.
In Fig. 7, we show the renormalization coefficient for the lowest mode as a function of the total number of modes that are coupled to the Josephson junction.Importantly, the renormalization coefficient appears to converge to a value slightly below 0.90.This supports our choice of the mode cutoff frequency, based on the superconducting gap, which yielded a renormalization coefficient of approximately ẼJ /E J ≈ 0.90 (see Sec. III in the main text).
To understand this behavior, we consider the renormalization coefficient expressed in the harmonic-oscillator basis as where λ m denotes the zero-point fluctuations of the mode m.In the exponent, both the inverse of the mode angular eigenfrequency and the effective charging energy are present, which are both depicted in the inset of Fig. 7 as functions of mode number m. Due to the CPW structure, the eigenfrequency scales as ∼ 1/m.This scaling alone would lead to E * J /E J ∼ 1/M , causing the renormalization coefficient to vanish in the limit M → ∞.However, the effective charging energy appears to scale as ∼ 1/m γ , where γ > 1.Note that there is a clear change in the scaling behavior around m ≈ 20, after which the effective charging energy decreases faster than the inverse eigenfrequency, implying that γ > 1 in the limit m → ∞.In the context of the exponent in Eq. (D1), this results in a scaling given by λ m ∼ (1/m) 1+γ .Should the assumption γ > 1 hold true, it implies that the renormalization coefficient converges to a nonzero value, consistent with the behavior observed in Fig. 7.
Similar arguments regarding the convergence of the renormalization effects, especially when considering a Josephson junction that is capacitively coupled to the end of a CPW, are discussed in Refs.[38,50,51].These works demonstrated that the magnitude of the coupling between the qubit mode and the high-frequency modes exhibits a natural cutoff frequency.This cutoff is exclusively dependent on the Josephson capacitance C J as M → ∞.We observed similar behavior in the mode-dependent effective charging energies of the unimon circuit.Specifically, in the inset of Fig. 7, a natural frequency cutoff is apparent around M ≈ 50.Although not shown in the figure, this cutoff is strongly influenced by the selected value of C J .

Appendix E: Accuracy of the multimode model
In this appendix, we address the accuracy of the multimode model compared with the auxiliary-mode model which are described by the Hamiltonians given in Eqs. ( 12) and (1), respectively.Figure 8 illustrates that significant discrepancies in the lowest-mode anharmonicity between the two models arise primarily at two flux bias points: Φ diff /Φ 0 = 0.39 and 0.50.We observe that if the number of auxiliary modes increases from M = 2 to M = 4, the anharmonicity undergoes corrections of approximately 15 MHz in magnitude.In both instances, these corrections move the value closer to what is obtained with the multimode model.In addition, the inset of Fig. 8 reveals a convergent behavior for the multimode model, exhibiting markedly smaller deviations than the auxiliary-mode model as the energy cutoff E cutoff surpasses 100 GHz.For instance, if E cutoff is increased from 150 GHz to 200 GHz, the corrections to the anharmonicity is roughly 1 MHz.
Appendix F: Dependence of the unimon physics on E J and Zc We investigate the interplay between the Josephson energy, E J , and the characteristic impedance, Z c , in a unimon operated at its sweet spot of φ = π.We aim to understand how changes in these parameters influence the behavior of the system.
Starting with the auxiliary-mode Hamiltonian given by Eq. ( 1 where we have used effective inductance defined in Eq. (A23).We postulate that, given the high plasma frequency ω p = 1/ √ L J C J compared to the frequencies of the primary normal modes, the Josephson capacitance C J in the effective capacitance can be reasonably neglected.This assumption is reinforced by our numerical findings.
The characteristic impedance Z c manifests in the relevant quantities as Applying these relations to Eq. (F1), it follows that the cosine term is the only one depending on Z c .An expansion using a Taylor series provides where β 2 and β 4 are constants, and φ = 2π Ψ/Φ 0 .This expansion underscores an essential observation, that E J and Z c influence the system in roughly identical ways up to the second order.However, the distinction emerges in higher-order terms.Notably, the sensitivity of anharmonicity to changes in Z c increases with increasing Z c .Yet, one must also account for the renormalization effects from other modes.These effects, which become more pronounced at greater Z c , act to temper the increase in anharmonicity.
In summary, our investigation reveals that, within the bounds of our assumptions, an equal relative change in either E J or Z c produces an identical outcome on the system up to linear order.The differences primarily arise in high-order behavior and in the distinct effects of renormalization at an elevated Z c .
is the half difference of the external magnetic fluxes shown in Fig.1(a), ϕ 0 is the dc magnetic flux offset across the Josephson junction, E J is the Josephson energy, Φ 0 denotes the flux quantum, Ω m is the resonance angular frequency of the auxiliary mode m, and {ξ m } are the coupling strengths of the auxiliary modes to the nonlinear mode.The definition for the effective lumped-element capacitance is

FIG. 3 .
FIG. 3. Energy levels of the unimon circuit as functions of (a) the location of the Josephson junction xJ, (b) Josephson energy EJ, and (c) characteristic impedance Zc of the CPW.The solid lines represent the qubit subspace of the three first modes, for which the rest of the modes are in the vacuum state.All other states are represented by dashed lines.The inset of panel (a) provides a close view of an example avoided crossing.For these simulations, we used the parameters from TableI, with the number of lower modes set to M0 = 8, except that the sweep of the characteristic impedance is conducted such that the frequency of the lowest auxiliary mode, given by Ω1 = π/(2l √ L l C l ), remains constant.
FIG. 4. (a) Anharmonicity of the unimon qubit and (b) the renormalization coefficient ẼJ/EJ as function of the location of the junction xJ.We show in (a) the anharmonicity obtained using the single-mode approximation (solid line), renormalization model (dashed line), and the multimode model (dotted line) detailed in Sec.II.The anharmonicities computed directly from the auxiliary-mode Hamiltonian at xJ = 0 are also indicated by markers, accounting for different numbers of auxiliary modes.The inset in (b) shows the contribution from each mode to the renormalization coefficient for xJ/l = 0 (green line) and xJ/l = 0.46 (orange line) [see Eq.(11)].These numerical calculations are carried out using the parameters listed in TableI.

FIG. 5 .
FIG. 5. Anharmonicities and cross-Kerr interaction energies for the two lowest modes of the unimon circuit as functions of (a) the location of Josephson junction xJ, (b) Josephson energy EJ, and (c) characteristic impedance Zc of the CPW.The dashed vertical line in panel (a) indicates the value of xJ, which is subsequently used for the graphs in panels (b) and (c).The notation αm(i, j) corresponds to the anharmonicity of the m:th mode (m = 1, 2, 3), with the other two modes having i and j excitations.The left argument (i) pertains to the lower mode of the two arguments.The cross-Kerr interaction is denoted by Knm(i), where n and m (n, m = 1, 2, 3) represent the modes between which the interaction is computed.The argument (i) equals to the number of excitations in the mode not primarily involved in the interaction.The simulation utilize the parameter values of TableI.

FIG. 6 .
FIG. 6. (a)-(c) Anharmonicities and (d)-(f) cross-Kerr interaction energies for the three lowest modes of the unimon circuit as functions of (a) and (d) the Josephson junction location xJ, (b) and (e) Josephson energy EJ, and (c) and (f) characteristic impedance of the CPW Zc.The notation αm(i, j) and Knm(i) is defined in Eqs.(15) and (16), respectively, and in Fig.5.For the the parameter values that are not swept, we use TableI, except for the Josephson junction location set at xJ/l = 0.62, as marked by the vertical dashed line in panels (a) and (d).

FIG. 7 .
FIG. 7.Renormalization coefficient ẼJ/EJ (M0 = 1) as a function of the total number of modes M (black curve).The inset shows the inverse of the mode eigenfrequency (red curve) and effective charging energy EC,m = c 2 m e 2 /(2C eff ) (blue curve) as functions of the mode number m.The calculations are carried out for xJ/l = 0 and ϕ0 = π.Note that only modes that couple to the Josephson junction are included.In this case only every other mode is coupled.

FIG. 8 .
FIG. 8. Anharmonicity of the lowest mode as a function of the Φ diff /Φ0 for the multimode (M0 = 10) and auxiliarymode (M = 2) models.In addition, two data points for the auxiliary-mode model with M = 4 are shown at positions where the deviation between the two models is most pronounced (Φ diff /Φ0 = 0.39 and 0.50).The inset shows the anharmonicity of lowest mode as a function of the energy cutoff E cutoff , as determined by the multimode model.The parameters used for these calculations are xJ/l = 0, L l = 0.83 µH/m, C l = 83.0pF/m, and EJ/h = 19.0GHz.
), we define new rescaled flux operators as Ψ = √ C eff Ψ and χm = √ C eff χm .The corresponding charge operators are defined as Q = Q/ √ C eff and Ξm = Ξm / √ C eff .This modification presents us with the

TABLE I .
Physical parameters used in the simulations of this paper unless otherwise explicitly stated.