Laser-annealing Josephson junctions for yielding scaled-up superconducting quantum processors

As superconducting quantum circuits scale to larger sizes, the problem of frequency crowding proves a formidable task. Here we present a solution for this problem in fixed-frequency qubit architectures. By systematically adjusting qubit frequencies post-fabrication, we show a nearly ten-fold improvement in the precision of setting qubit frequencies. To assess scalability, we identify the types of 'frequency collisions' that will impair a transmon qubit and cross-resonance gate architecture. Using statistical modeling, we compute the probability of evading all such conditions, as a function of qubit frequency precision. We find that without post-fabrication tuning, the probability of finding a workable lattice quickly approaches 0. However with the demonstrated precisions it is possible to find collision-free lattices with favorable yield. These techniques and models are currently employed in available quantum systems and will be indispensable as systems continue to scale to larger sizes.


I. INTRODUCTION
Realizing robust large-scale quantum information processors is one of the foremost challenges in quantum science. Many practical applications have been proposed for robust quantum computers, including estimating the ground state energy of chemical compounds and implementing machine learning algorithms [1][2][3][4][5]. Quantum advantage relative to classical computers can be realized without full fault-tolerance, but requires large quantum circuits that a classical computer cannot simulate [6]. Recent demonstrations have shown qubit circuits nearly at the threshold for demonstrating quantum advantage [7]. Much work remains in order to realize fault-tolerant quantum processors; however, scale-up of solid-state quantum circuits has shown consistent and ongoing progress [8][9][10][11][12][13][14][15]. As the qubit circuits are scaled up, they must maintain high oneand two-qubit gate fidelities, high qubit connectivity, and low cross-talk error which can be measured in a holistic sense via the quantum volume of the circuit [16]. Lattices of fixed-frequency transmon qubits represent a promising architecture for building systems of larger sizes [7]. A growing number of systems at the 20 to 50-qubit scale are now available to users through cloud access. A variety of technical challenges confront further system scaling, including improving 3-dimensional circuit integration and qubit coherence. High on the list of such challenges is the issue of 'frequency crowding.' Fixed frequency transmon qubits using the two-qubit cross-resonance (CR) gate form a promising architecture for scaling up quantum systems. Fixed-frequency transmons are largely insensitive to charge or flux noise, and have achieved coherence times of 100 µs and growing. The CR gate, a hardware-efficient all-microwave gate [17][18][19][20], is readily used to entangle these qubits with gate fidelities above 99%, approaching the threshold for fault-tolerant codes [21]. To achieve these fidelities, the CR gate needs not only high coherence qubits, but also a precise setting of the qubits' frequencies. The CR gate activates a ZX interaction by driving one 'control' qubit with a microwave pulse at the other 'target' qubit's transition frequency. The magnitude of the ZX as well as other Hamiltonian terms depends on the relative frequencies of the two qubits [22,23]. Diminished ZX magnitude increases gate time, while other terms such as ZZ add gate errors. Neighboring qubits having the wrong detuning will exhibit a 'frequency collision' in which the ZX may be suppressed or other undesirable effects arise.
Maintaining high gate fidelities for all pairs in a lattice will require solving this 'frequency crowding' problem by precise setting of qubit frequencies to specified values, as characterized by a standard deviation σ f . To achieve low σ f , the tunnel-junction conductance must be controlled with high precision. Transmon frequency f 01 follows hf 01 √ 8E J E C − E C , where Josephson energy E J = Ic 2e is many times greater than charging energy E C = e 2 2C [24]. In typical transmons, a photolithographically defined capacitance C has dimensions in the tens to hundreds of microns and varies little from qubit to qubit. The critical current I c is set by a tunnel barrier of area ∼ 100 × 100 nm and thickness a few nm, and is thus challenging to fabricate with precision better than a few percent [25][26][27][28]. However, tunnel barrier resistance R n is readily measurable to precision better than 0.1% and relates to I c according to the Ambegaokar-Baratoff relation I c = π∆ 2eRn (where ∆ is the superconducting gap energy) [29]. We can therefore measure R n before a chip is cooled in order to assess qubit frequency imprecision. The best demonstrated precision in setting R n at time of fabrication is 2% [28]. A 2% variation in R n indicates a fractional σ f of 1%.
Careful design of lattices can enable error correction codes while at the same time minimizing the likelihood of 'frequency collisions' and therefore the required σ f for fabrication yield [30,31]. Yet even the most robust designs require a fractional σ f of 0.25% to 0.5%, which represents a factor of 2 to 4 improvement over the 2 best literature results. To overcome such limits will require rework of individual qubits' tunnel junctions after fabrication. Thermal anneal has been shown to increase tunnel resistance R n , and laser heating has been demonstrated as a highly localized re-work tool [32][33][34][35][36][37]. However, the inherent variability of the anneal process itself must be overcome, and qubit frequency control utilizing such techniques at scale has never been presented in the literature.
In this paper, we introduce an adaptive post-fabrication trimming technique that we use to incrementally adjust R n on a qubit-by-qubit basis, thereby overcoming inherent variability in both initial qubit fabrication and the laser anneal. For the first time, an improvement in qubit frequency precision is demonstrated in terms of narrowed frequency distributions. Crucially, we demonstrate qubit frequency imprecision σ f of the same magnitude as the imprecision of predicting f 01 from R n . To estimate the scalability of this technique for the fabrication of error-corrected lattices, we employ a statistical yield model based on σ f relative to specific collision bounds. This model predicts the severity of the frequency crowding problem for different topologies and scales of error corrected multi-qubit lattices as a function of code distance. The model demonstrates that using conventional transmon fabrication, scaled-up qubit lattices will fail to evade 'frequency collisions'. However, our novel trimming technique achieves adequate σ f for scalable fabrication of distance-3 through distance-7 heavy-square and heavyhexagon codes. In particular, this technique enables the high yield fabrication of the distance-3 and distance-5 heavy-hexagon lattices currently deployed as IBM cloud connected systems [38].

A. Frequency Precision σ f From Transmon Fabrication
To assess the σ f resulting from qubit fabrication, we developed a test vehicle containing a large number of identically-fabricated qubits ( fig. 1). We cooled the chip in a dilution refrigerator and used dispersive readout through half-wave microwave resonators to measure qubit frequencies [39]. We measured the frequencies of 31 qubits to a precision better than 100 kHz using a Ramsey fringe method. The qubit frequencies had random variation σ f = 132.3 MHz ( fig. 1), or 2.3% of the median frequency. After warming the qubits to room temperature, we measured their junction resistances. In fig. 1, we show a plot of R n compared to transmon frequency, demonstrating that the observed variation in σ f is accounted for almost entirely by R n variation. The behavior may be fit to a power law of approximately − 1 2 power, as expected from transmon theory and the Ambegaokar-Baratoff formula.
For a population of transmon qubits whose frequency scatter is dominated by scatter in R n , we expect the fractional standard deviation in R n to be twice that of σ f . This is consistent with the standard deviation in junction resistances which is found to be σ R of 365 Ω, or 4.6% of the median R n . We also assess the fidelity of the frequencies to the f -vs-R n correlation in terms of the residual scatter after subtracting the fit line from the frequency values. This appears in the inset in fig. 1 and exhibits a standard deviation 14.5 MHz, or 0.25% of the qubit median frequency.

B. Tuning Using Selective Laser Anneal
To reduce σ f , we developed a technique for selective laser anneal to shift tunnel resistance R n by pre-calibrated increments. (See section IV and fig. 2). We demonstrate the achievable frequency control of this technique by shifting the 31 measured qubits into a two-frequency pattern. We employed an R n vs f correlation ( fig. 1) to designate the target resistances. We shifted 16 junctions to one target R n and 15 to another target R n . After tuning, the group of 16 junctions had median resistance 7.984 kΩ and the group of 15 had median resistance 8.798 kΩ. The 31 junctions reached their targets with an overall precision of σ R = 51Ω, about 0.61%. In a dilution refrigerator, we re-measured the frequencies of the qubits in the two groups. Aside from two of the qubits, which we measured using CW spectroscopy (precision 2 MHz), all qubits were remeasured in the same way as in the first cooldown. The resulting frequencies appear in fig. 1. The two frequency groups are approximately normally-distributed and have medians f 0,1 = 5.430 GHz and f 0,2 = 5.7046 GHz. Cal- where f 0,j represents f 0,1 or f 0,2 as appropriate for a given qubit Q i , we assess the overall precision σ f = 14.0 MHz. This imprecision is nearly identical to the residual scatter from the f (R) fit line ( fig. 1) which guided the tuning, and the fractional precision σ f / f = 0.25% is slightly better than half of the fractional precision in setting R n . By these comparisons, we see that in this experiment σ f is limited by both the precision of setting R n and the precision of predicting f from R n . Drift in R n reported in the literature [32] does not appear to be a limiting factor in this study. Achieving smaller σ f will require improvements in setting R n . As we show in section IV, the laser-anneal tuning technique is capable of precisions of 0.3% in R n . On that basis the imprecision of 14.5 MHz in predicting f from R n would dominate the imprecision in σ f .

III. DISCUSSION
Our post-fabrication trimming reduced σ f by 9.5× compared to initial fabrication. To assess whether this level of precision is sufficient to reliably prepare lattices of fixed-frequency transmons capable of error-correcting codes, we must quantify the frequency-crowding problem. Transmon qubits are weakly anharmonic and have decreasing transition energies at higher levels. Therefore, degeneracies among the |0 → |1 , |1 → |2 and |0 → |2 transitions of nearby qubits can all contribute to 'frequency collisions. ' We must consider the relative frequencies of both nearest-neighbors and next-nearest-neighbors in the lattice [22,41,42]. Fig  3 illustrates the relative positions of nearest-neighbor and next-nearest-neighbor qubits in a section of lattice, and table I lists the seven cases most likely to lead to gate errors [22]. We can think of them qualitatively as follows: Type 1 causes hybridization of states in Q j and Q k , while in type 2 the CR pulse excites Q j into the non-computational |2 state. Type 3 excites Q k to the |2 state, but does not require a CR tone. In condition 4, ZX is weak, which implies long gate times and increased gate error [22,23]. In type 5, the CR gate addresses an additional neighboring qubit. In type 6, when one qubit is the target of a CR gate, its next-nearest-neighbor leaks to the |2 state. Type 7 causes Q j to leak to the |2 state during a CR gate.
Around each of the 'frequency collisions' described in Existing multi-qubit systems with CR gates typically exhibit two-qubit gate errors of 1 to 2 % regardless of frequency [12,41]. Ref [22] considers an effective-Hamiltonian model for the CR gate, as a function of the relative frequency of control and target qubits. We use this model to estimate the frequency windows for nearest-neighbor collisions (table I, types 1 to 3) to cause gate errors exceeding ∼ 1 %. We make an assumption that similar bounds apply to next-nearest-neighbor interactions (types 5 to 7).
A useful lattice of qubits should enable high quantum volume and fault-tolerant operation while avoiding all of the 'frequency collisions' and forbidden regions presented in table I. Both lattice layout and the pattern of qubit frequencies are relevant. We consider three types of lattices: square, 'heavy-square' and 'heavy-hexagon' (fig. 3). Lattices comprise qubits and two-qubit connections, each qubit being linked to no more than four neighbors. In many practical implementations, these links comprise microwave-resonant buses. A square lattice facilitates 'surface code' fault-tolerant codes [43]. Recent literature describes hybrids of the surface code with Bacon-Shor type codes, which can be employed in 'heavy hexagon' and 'heavy square' lattices to achieve fault-tolerance, albeit with lower error thresholds than the surface code [30]. In addition to the data and ancilla qubit roles employed in the surface code, these hybrid codes assign a portion of the lattice as 'flag' qubits.
In the square lattice, every qubit in the bulk of the lattice lies on a degree-four vertex, while some at edges have degree two or degree three. If we populate the

Type Definition
Participants Nearest-neighbor qubits Q j , Q k ± 30 MHz 4 f k,01 < f j,12 or f j,01 < f k,01 Control qubit Q j , target qubit Q k -5 f i,01 = f k,01 Q j is control to Q i and/or Q k & is nearest-neighbor to both. ± 17 MHz 6 f i,01 = f k,12 or f i,12 = f k,01 Q j is control to Q i and/or Q k & is nearest-neighbor to both. ± 25 MHz 7 f j,02 = f k,01 + f i,01 Q j is control to Q i and/or Q k & is nearest-neighbor to both. ± 17 MHz  fig. 3. Bounds for types 1,2,3 are estimated from model results [22] as the region where gate errors due to 'frequency collision' exceed ∼ 1 %. Bounds for types 5,6,7 are based on those in types 1 and 3.
square lattice with 5 distinct frequencies of qubits, In contrast to the square lattice, the 'heavy square' lattice includes both degree-two and degree-four vertices in the bulk. Degree-one, -two or -three vertices appear at the edges. We take advantage of this pattern to make all the degree-two vertices control qubits, using a three-frequency pattern f 3 > f 2 > f 1 . Since every control qubit (frequency f 3 ) is linked to at most two target qubits, we need only two properly-chosen target-qubit frequencies (f 1 and f 2 ) to satisfy conditions 5, 6 and 7 of table I, as shown in fig. 3. A third type of lattice, the 'heavy hexagon', uses a similar scheme. Here the bulk of the lattice includes degree-three and degree-two vertices. Additional degree-two and degree-one vertices lie at the edges. In this lattice, all of the 'frequency collisions' and forbidden regions can be satisfied using only three frequencies f 3 > f 2 > f 1 , with all control qubits residing on degree-two vertices with frequency f 3 .
collisions'. As σ f increases, the number of 'frequency collisions' rises steadily. As σ f → f 01 − f 12 , the different conditions appearing in table I all become likely, and a limiting number of 'frequency collisions' is reached. Yield follows the inverse trend, as seen in fig. 4. As σ f increases, the likelihood of finding a 'collision free' chip falls off sharply. While the step sizes between frequencies f 1 to f 5 are important, absolute frequency values are not. Setting f 1 = 5.0, f 2 = 5.07 and f 3 = 5.14 GHz works as well as f 1 = 5.05, f 2 = 5.12 and f 3 = 5.19 GHz.
The yield and mean collision number are a function of the several different collision types and bounds, so they are not readily susceptible to an analytic formulation. However, we can propose a simplified model for yield: in order for a lattice to be collision-free, every qubit in the lattice must fall within some frequency 'window' ±∆f relative to its setpoint. Presuming the qubit frequencies are normally distributed, the probablity of this occurring goes as the cumulative distribution function, raised to the power N , where N is the number of qubits: In the yield plot in fig. 4, we fit this expression to find ∆f for each lattice.
These model results allow us to predict how different lattice types and frequency patterns will respond to fabrication imprecision. As shown in fig. 4, if imprecision σ f is greater than 30 MHz, any d = 5 lattice will exhibit > 10 'frequency collisions' of one or another of the types listed in table I, causing the affected gates to have error rates above ∼ 1 %. However, if σ f = 10 MHz then on average the d = 5 square lattice will exhibit 5 'frequency collisions', while the 'heavy square' and 'heavy hexagon' lattices will exhibit 0.1 'frequency collision'. Considered in terms of yield, we see from fig. 4 that if σ f = 10 MHz, then for a d = 5 device, a square lattice with 5-frequency pattern has a 0.8% likelihood to be 'collision free', whereas a 'heavy square' lattice with 3-frequency pattern has 90% likelihood and 'heavy hexagon' with 3-frequency pattern has 92% likelihood. Alternatively we can ask, how well do we have to control σ f ? If we seek a 10% yield, then fig. 4 indicates that for a d = 5 device, a square lattice with 5-frequency pattern requires σ f < 8 MHz, whereas a 'heavy square' lattice with 3-frequency pattern requires σ f = 16 MHz and 'heavy hexagon' with 3-frequency pattern requires σ f = 17 MHz. Although the square lattice requires 10 to 20 % fewer qubits than the other types at each distance d, it requires far better frequency precision.
The as-fabricated σ f seen in fig. 1  As seen from the Monte Carlo analysis, the laser-anneal rework method can scale to the > 100 qubit size, enabling a well-chosen lattice and frequency pattern to implement d = 7 error-correction codes free of frequency-crowding. To examine needs for the next generation of chips up to the 1000-qubit level, we can coarsely estimate requirements by extrapolating the fixed window model for the  fig. 3, capable of distance-5 codes, as well as to d = 3 and d = 7 scale lattices of square, heavy-square or heavy-hexagon type. (See Supplementary Figures 1-9 for lattice layouts and table II for numbers of qubits.) Color-coded dotted lines in part (b) are fits of each lattice yield to expression where N is the number of qubits and ±∆f defines an allowable 'window' around frequency set-points. (See table II.) Using this expression, two solid red lines predict yield for the heavyhexagon lattice type at 300 and 1000 qubits, using ∆f = 27.99 and 26.32 MHz, respectively. See fig. 5 for estimation of ∆f as a function of qubit number. heavy-hexagon lattice as shown in fig. 5. While the σ f = 14.0 MHz demonstrated here enables practical yield up to the 100 to 200 qubit scale, it is clear that roughly a factor of two further improvement is needed to scale towards 1000 qubits. Since this precision is also better than the resistance-to-frequency prediction precision shown in this work, development of further refinements in tuning and frequency prediction approaches will be necessary as the scale of fixed-frequency transmon circuits surpass the 100 qubit milestone.

A. Chip Fabrication
A chip of the kind used to determine σ f and to test our laser-anneal rework process appears in fig. 1. All microwave elements comprise Nb films ∼ 200 nm thick on a silicon substrate. Each qubit is coupled to a readout resonator but is not directly coupled to any nearby qubits. All transmon capacitors are identical. Junctions are fabricated using identical electron-beam lithographic patterns and deposited simultaneously using double-angle deposition and oxidation [44]. The individual qubit design is similar to that used in Ref [21] with aharmonicity f 12 − f 01 -330 MHz. Junctions have linear dimension ∼ 100 nm and are designed for I c of ∼ 30 nA. During packaging, we accidentally damaged three of the 36 qubits and found these to be non-functional when cooled in a dilution refrigerator. We left two of the remaining 33 qubits un-tuned as experimental controls, so that our tuning demonstration includes 31 qubits.

B. Tuning Using Selective Laser Anneal
We have built an integrated junction rework system that can measure and modify the junction resistance. Fig. 2 shows a schematic of our laser annealing system, which we call "Laser Annealing of Stochastically Impaired Qubits" (LASIQ). The laser output is generated by a diode-pumped solid-state laser, frequency doubled to 532 nm. Active power control of the anneal beam is performed using a piezo-rotary mounted waveplate and polarizing beam splitter (PBS), which is adaptively adjusted based on a pick-off beam measured on a downstream silicon photodiode. A precision-timed shutter is  . 1) in different lattices and frequency patterns ( fig. 3 and Supplementary Figures 1-9). The ∆f values correspond to the fit lines in fig. 4: in order for the device to be collision-free, every qubit in the lattice must fall within ± this 'window', relative to its set-point.
used to control the anneal duration, and beam alignment is performed using a mechanical mirror mount which directs the beam via pattern recognition to the transmon junction center. The beam is shaped as needed to avoid illuminating the junction directly [40].
By careful control of laser power and pulse duration, we use this system to adjust R n . This process overcomes the imprecision due to transmon fabrication, with a residual imprecision σ f due to the rework process. To develop the process, we prepared a set of more than 150 junctions identically to qubit junctions, and measured their response to a range of laser powers and exposure times. We recorded R n shifts up to 15% relative to initial R n , for anneal durations varying by an order of magnitude and laser powers varying by 20%. Response to laser power in particular was highly nonlinear. Based on these empirical calibrations of R n shift to power and exposure, we established a qubit tuning process: We first measure the transmon junction's R n using four-point probing of the transmon capacitor pads at 25 • C. Using a f (R n ) prediction based on a previously determined correlation curve ( fig. 1), we assign the junction a target resistance corresponding to the target frequency in a multiqubit chip lattice. Because the anneal can shift R n in only one direction, the target must be higher than the initial R n . We anneal the qubit junction using laser power and duration chosen from our calibration set, then re-measure its R n . A junction requiring large shifts in R n may require repeated anneals to reach its target, as shown in fig. 2. The control algorithm increases the resistance until the measured value is within 0.3% of the target value. In a separate trial of tuning precision, more than 300 junctions were tuned to target R n s ranging from 0.4% to 14.5% above their initial values, and landed successfully within this 0.3% margin. We expect 0.3% imprecision in R n to introduce 0.15% imprecision in transmon frequency.

C. Monte Carlo Frequency-Crowding Model
Using a Monte Carlo model, we can estimate the incidence of 'frequency collisions' in a lattice as a function of σ f . We assume that imperfect frequency-setting will distribute qubit frequencies normally around their design frequencies with standard deviation σ f . For lattices of the type shown in fig. 3, we designate 3 to 5 frequencies f 1 , f 2 , f 3 , f 4 , f 5 spaced at regular intervals in the pattern shown. We set f 1 = 5 GHz, similar to real-world transmons [38,45]. We sample the qubit frequencies randomly around these values and count the collisions throughout the lattice, as listed in table I. This process is illustrated in fig. 6. We repeat the frequency-assignment and counting to build statistics for a given lattice and frequency pattern. We then repeat the model for a range of σ f values from 0 to 150 MHz. We repeat the entire process over a range of frequency spacings, to find the spacing that minimizes 'frequency collisions' at each value of σ f . As a function of σ f we can then extract 1) the mean number of total collisions in the lattice, and 2) the fraction of repetitions which result in zero collisions ('yield'). Our simulations used 1000 repetitions except to find yield below 1% in d = 5 lattices and below 0.2% in d = 3 lattices, which used 4000 repetitions, and in d = 7 lattices to find mean collisions for σ f < 16 MHz or yield above 50% (100 repetitions) or to find mean collisions for σ f > 16 MHz (40 repetitions fig. 3: f1 = 5.00, f2 = 5.07, f3 = 5.14, f4 = 5.21, f5 = 5.28 GHz. To model the lattice statistically, treat the frequencies f1 to f5 as means of distributions. Right: Mean frequencies and normal distributions characterized by σ f . For each position in the lattice, sample from the local distribution. Choose random frequencies in this fashion, count collisions as described in table I and repeat to gather statistical sample. Process is repeated for differing spacings between mean frequencies f1 to f5, and for different distribution widths σ f .