Pulse optimization for high-precision motional-mode characterization in trapped-ion quantum computers

With the advent of promises to deliver a quantum computer of scale beyond the capability of conventional computers, comes the responsibility to indeed manufacture one that can reliably be tuned for its normal operations. Without the ability to tune, a quantum computer may produce unpredictable and erroneous results. This necessitates the development of methods that can efficiently and accurately characterize the system parameters, as a prerequisite for its high-fidelity operations.

In this paper, we aim to address this challenge for a trapped-ion quantum computer. In particular, our focus is to enhance the accuracy of motional-mode parameter estimation Kang et al. (2023a), as these parameters are essential for designing and performing high-fidelity two-qubit operations Cirac and Zoller (1995); Mølmer and Sørensen (1999); Sørensen and Mølmer (1999); Blümel et al. (2021, 2021); Li et al. (2022). An accurate characterization of mode parameters can provide several additional benefits as well. First, the overhead in gate calibrations caused by incorrect parameter estimates Maksymov et al. (2021); Gerster et al. (2022) can be reduced. Second, gate pulses can be more efficient in terms of control-signal power, gate duration, and robustness than pulses designed with inaccurate motional-mode parameters Blümel et al. (2021, 2021).

In a characterization experiment, a model is used to predict the physical parameters being characterized via fitting to some measured quantities. A major bottleneck to achieving high-accuracy characterization is the mismatch between the model’s predictions and experimental outcomes Kang et al. (2023a), resulting from various sources of experimental imperfections and/or model violations. Ref. Kang et al. (2023a) presents improved fitting models as well as parallelized characterization protocols to address some of these mismatches. To achieve further improvement in the accuracy of mode-parameter characterization, here, we investigate a complementary approach to the characterization tools developed in Kang et al. (2023a) by leveraging the control degrees of freedom that are already readily available for the quantum gate implementation on the hardware level.

Specifically, we use pulse optimization, which is a widely used tool in quantum optimal control Werschnik and Gross (2007). Indeed, many pulse-shaping paradigms have been developed, in particular for the purpose of designing high-fidelity single-qubit and two-qubit gates Khaneja et al. (2005); Caneva et al. (2011); Maximov et al. (2008). For trapped ions, the pulse shaping techniques are standard and have been previously adopted to accomplish amplitude-modulated (AM) Zhu et al. (2006); Roos (2008); Kim et al. (2009); Choi et al. (2014); Debnath et al. (2016); Figgatt et al. (2019), frequency-modulated (FM) Leung et al. (2018); Landsman et al. (2019); Wang et al. (2020); Kang et al. (2021, 2023b), and phase-modulated (PM) Shapira et al. (2018); Milne et al. (2020); Green and Biercuk (2015); Lu et al. (2019) entangling gates. Here, we devise pulse shaping paradigms to eliminate the effect of dominant error sources in the characterization process, which in turn can help improve the performance of multi-qubit gates Katz et al. (2023a, b) and parallel gates Figgatt et al. (2019); Bentley et al. (2020); Grzesiak et al. (2020). We note that a recent work explored pulse shaping for system parameter characterization in a single-qubit setup Stace et al. (2022).

The rest of the paper is structured as follows. In Sec. II, we provide necessary preliminary information for concretely defining our mode characterization problem. We then formalize the problem in Sec. III, where we focus on the cross-mode coupling and detuning errors, to be detailed in the section, as the two dominant sources of error in the mode characterization. In Sec. IV, we explore pulse shaping, where we exploit the degrees of freedom in modulating the beams that illuminate the ions to engineer the response of the system, such that the previously discussed two errors are minimized. In Sec. V we demonstrate the viability of our proposed pulse shaping schemes via extensive numerical simulations on a three-ion chain. Finally, we discuss the applications and potential limitations of our proposed scheme in Sec. VI and provide outlook in Sec. VII.

In a typical trapped-ion quantum computer, a qubit is encoded in two internal states of an ion. As $N$ ions form a linear Coulomb crystal, the collective motion of the ions can be decomposed into $3N$ normal modes, each of which is described as a quantum harmonic oscillator. A pair of qubits within the ion chain can be entangled via a state-dependent force, frequently using lasers. Such force couples the ions’ internal states to the motional modes, typically strongly coupling to a subset of $N^{\prime}$ modes, where the motional modes serve as a “bus” that establishes the communication between the target pair of ions that we want to entangle.

To produce an ideal entangling gate, two criteria must be satisfied; first, the residual entanglement between the qubits and the motional modes must be removed at the end of the gate; second, the degree of entanglement between the target pair of qubit states must reach a pre-specified amount. Indeed, accurate knowledge of the motional-mode parameters, mode frequencies and mode-ion coupling strengths (Lamb-Dicke parameters), is necessary to exactly satisfy the criteria. In practice though, the parameters are inevitably known only to finite accuracy, due, e.g., to faulty characterization or experimental drifts. This necessitates gate pulse design methods to construct solutions that are robust to parameter offsets, albeit at the cost of a significantly higher power or longer pulse length Blümel et al. (2021, 2021); Kang et al. (2023b); Shapira et al. (2018); Milne et al. (2020); Jia et al. (2023); Valahu et al. (2022). Despite much effort, readily available rough theoretical estimates of the mode parameters, calculated using minimal knowledge of the system configurations such as the voltages applied to the trap Wineland et al. (1998), have been proven to be insufficient for high-fidelity quantum operations, as evidenced in Ref. Blümel et al. (2021) (see also footnote [35] of Ref. Kang et al. (2023a)).

It is thus not surprising that spectroscopic methods Mavadia et al. (2014); Goodwin et al. (2016); Stutter et al. (2018); Hrmo (2018); Welzel et al. (2018); Joshi et al. (2019); Hrmo et al. (2019); Jarlaud et al. (2020); Feng et al. (2020); Chen et al. (2020); Sosnova et al. (2021) have been employed to better estimate the mode parameters. Briefly, in all these methods, one starts by initializing the ion’s internal and motional state to the ground state $|0,0\rangle_{j,p^{*}}$ following the standard initialization and cooling procedures. Here $|a,b\rangle_{j,p^{*}}$ denotes the composite state of the computational basis state $|a\rangle_{j}$ of ion $j$ with $a\in\{0,1\}$ and the motional Fock state $|b\rangle_{p^{*}}$ of mode $p^{*}$ with phonon number $b$, and $p^{*}$ is the index of the mode we aim to characterize. One then turns on the Hamiltonian

$\displaystyle\hat{H}_{I,j}(t)$

$\displaystyle=\hat{\sigma}^{+}_{j}\exp{i\sum_{p=1}^{N^{\prime}}\eta_{j,p}(\hat{a}_{p}e^{-i\omega_{p}t}+\hat{a}_{p}^{\dagger}e^{i\omega_{p}t})}g_{j}(t)+h.c.,$

(1)

by illuminating ion $j$ for a duration $\tau$ with the pulse form $g_{j}(t)$ that in principle can be of any shape, where $\hat{\sigma}^{+}_{j}:=\left|1\right>\left<0\right|$ acts on the computational basis state of ion $j$, $\hat{a}_{p}^{\dagger}$ and $\hat{a}_{p}$ are the creation and annihilation operators acting on mode $p$, $\omega_{p}$ is the frequency of mode $p$, and $\eta_{j,p}$ is the Lamb-Dicke parameter that couples ion $j$ and mode $p$, defined according to $\eta_{j,p}:=b_{j,p}|\vec{k}_{j,p}|/\sqrt{2m\omega_{p}},$ where $b_{j,p}$ is the $j$-th element of the normalized eigenvector of mode $p$, $m$ is the ion mass, and $\vec{k}_{j,p}$ is the wavevector of the electric field that couples ion $j$ to mode $p$, projected along the motional direction of mode $p$. Here we assume $|g(t)|$ to be much smaller than the qubit frequency splitting.

Indeed in Mavadia et al. (2014); Goodwin et al. (2016); Stutter et al. (2018); Hrmo (2018); Welzel et al. (2018); Joshi et al. (2019); Hrmo et al. (2019); Jarlaud et al. (2020); Feng et al. (2020); Chen et al. (2020); Sosnova et al. (2021) $g_{j}(t)$ was chosen to be near-resonant to one of the blue-sideband (BSB) transition frequencies $\omega_{p^{*}}$, i.e., $g_{j}(t)={A}_{j}\exp\{-i(\mu_{j}t+\phi_{j})\}$ with the detuned drive frequency $\mu_{j}=\tilde{\omega}_{j}-\omega_{j}^{\rm{qbt}}$ $\approx$ $\omega_{p^{*}}$, where $A_{j}$ is the Rabi frequency of the resonant qubit-state transition, $\phi_{j}$ is the laser phase, $\tilde{\omega}_{j}$ is the frequency of the laser’s electric field, and $\omega_{j}^{\rm{qbt}}$ is the frequency separation between the two internal states of qubit $j$. The ion-mode system then undergoes the so-called BSB transition, where a Rabi flop between the states $|0,0\rangle_{j,p^{*}}$ and $|1,1\rangle_{j,p^{*}}$ occurs. Since the final $|1\rangle_{j}$-state population of ion $j$ at time $\tau$ has dependence on $\omega_{p^{*}}$ and $\eta_{j,p^{*}}$ according to (1), one can then extract the values of the mode parameters by measuring the population through statistics accumulated over multiple shots and comparing the theoretical predictions implied by (1) with the measurements.

We consider the aforementioned, conventional spectroscopy method using $g_{j}(t)={A}_{j}\exp\{-i\mu_{j}t+i\phi_{j}\}$ as our baseline, with which we compare our characterization method using pulse shaping, to be developed in later sections. Note the baseline uses a single-tone pulse, where the drive frequency $\mu_{j}$ is near-resonant with the frequency $\omega_{p^{*}}$ of the target mode $p^{*}$. We denote such non-modulated pulses as square pulses.

In theory, if one can simulate (1) efficiently and exactly, the mode parameters $\omega_{p}$ and $\eta_{j,p}$ for all $j$ and $p$ can in principle be estimated to within the model violation. However, such a simulation becomes a computationally prohibitive task as the system size increases. Traditionally, the so-called single-mode BSB model Wineland et al. (1998) has thus been employed to enable a scalable simulation, effectively focusing on a single mode and a single ion only, at the cost of introducing further, non-insignificant model-violation errors.

In practice, additional considerations need to be given. Recall the agreement between the model predictions and the experimental results in the BSB transition population is at the core of the mode-parameter estimation and the population has dependencies on both $\omega_{p}$ and $\eta_{j,p}$. Inaccurate estimation of $\omega_{p}$, a known parameter to fluctuate and drift Kang et al. (2023b), would thus result in inherent inaccuracy in estimating $\eta_{j,p}$ and vice versa.

To resolve the theoretical and practical issues then, two complementary approaches may be considered. The first approach involves finding more sophisticated theoretical models that better approximate the system while still remaining computationally feasible and developing tailor-made experimental protocols that better separate the effect of the uncertainties in $\omega_{p}$ from that of the uncertainties in $\eta_{j,p}$ in the various BSB populations induced empirically. This approach has indeed been studied in detail in Ref. Kang et al. (2023a), however using a square pulse. In a contrasting and complementary approach, which is the one explored and investigated here, active “noise canceling” can be considered, via pulse shaping: By deliberately constructing a non-square, shaped pulse to eliminate significant model-violation terms or desensitize the BSB-transition population with respect to parameter change, significant improvement in the mode-parameter estimation can be achieved.

In this section, we concretely lay out the necessary technical details for understanding our pulse shaping method and its benefits. By showing the approximations used to arrive at the simple, single-mode BSB model, we can pinpoint the terms responsible for the model-violation errors. Further, we can explicitly express the sensitivity term. Note, in the forthcoming discussion, for simplicity, we consider the case of illuminating one ion at a time for one target mode, thus omitting the reference to ion $j$ wherever contextually clear; Parallelization by simultaneously illuminating multiple ions and targeting multiple modes Kang et al. (2023a) may be considered but is beyond the scope of this paper.

To start, we reiterate that, when the mode characterization is performed for a trapped-ion quantum computer, experimental evolution is compared to that implied by the simulated Hamiltonian. Specifically, a series of experimentally observed BSB-transition population $P$ is compared with that expected from the simulation, denoted ${\mathcal{P}}$, to reveal the mode parameters of interest. As a result, the discrepancy between $P$ and $\mathcal{P}$ directly leads to errors in characterizing the mode parameters.

Multi-mode and Single-mode BSB Hamiltonian – To obtain $\mathcal{P}$ at scale beyond a handful number of ions, we first apply a series of approximations to the Hamiltonian in (1) to arrive at the single-mode BSB Hamiltonian, better amenable to a scalable simulation. Assuming small $\eta_{p}({\leftarrow}\eta_{j,p})$, we linearize the Hamiltonian in (1) by performing a first-order Taylor expansion on the exponential term in (1). Then, a series of rotating-wave approximation (RWA) is applied to remove the frequency components that are far off-resonant from all mode frequencies (see Appendix B for details). Here we assume a more general form of the pulse function $g(t)$ that have frequency componenet(s) near the mode frequencies only. These approximations combine to give the $N^{\prime}$-mode BSB Hamiltonian

$$\hat{H}^{\prime}_{I}(t)=i\sum_{p=0}^{N^{\prime}-1}\eta_{p}e^{i\omega_{p}t}g(t)\hat{\sigma}^{+}\hat{a}^{\dagger}_{p}+h.c.$$

(2)

Next, as only mode $p^{*}$ is targeted, we may apply a stronger RWA to further remove the frequencies of the non-target modes, assuming $|\eta_{p}A|\ll|\omega_{p}-\mu|\;\forall p\neq p^{*}$. This results in the single-mode BSB Hamiltonian

$$\hat{H}^{\prime}_{I,p^{*}}(t)=i\eta_{p^{*}}e^{i\omega_{p^{*}}t}g(t)\hat{\sigma}^{+}\hat{a}^{\dagger}_{p^{*}}+h.c.$$

(3)

See detailed derivations and approximations used to arrive at (2) and (3) in Appendix B. We note that the simpler model in (3) demands far less computational cost in its simulation, whereas the cost for simulating the more complicated model in (2) grows exponentially with the number of modes $N^{\prime}$. Nonetheless, the second RWA that we apply leading to (3) no longer holds when $|\eta_{p}A|$’s are of comparable magnitudes to the finite mode-frequency spacings $|\omega_{p^{*}}-\omega_{p}|$, which constitutes the heart of the problem that we describe next.

We are now ready to discuss our problem details that base the comparison between $P$ and ${\mathcal{P}}$, both of which we obtain by simulation in this paper. Starting with $P$, recall the most complete theoretical model considered herein is (1). Comparing (1) and (2), one can see that the leading-order difference is $O(\eta^{2})$, which is indeed the well-known Debye-Waller (DW) effect Wineland and Itano (1979); Wineland et al. (1998). Briefly, it affects the population $P$ by reducing the effective BSB Rabi frequency – the spread of the ion’s position wavepacket, captured in (1) but not in (2), leads to this reduction. Since our goal is to silence the leading-order, significant model-violation terms, which we show to be $O(\eta)$ in the paragraphs to come later, for our purposes, it suffices to consider the evolution implied by (2) as our simulated experiment that generates the (simulated) experimental population $P$. Further, the DW effect is known to be re-capturable by using a more advanced model that predicts the BSB population without adding significant computational cost Kang et al. (2023a). As for the model-based, approximate BSB population ${\mathcal{P}}$, we obtain it from simulating the evolution implied by the Hamiltonian in (3) over duration $\tau$.

Elimination targets for pulse shaping – The source of discrepancy between $P$ and ${\mathcal{P}}$ forms the elimination target by our pulse shaping. Such an elimination is in general not achievable via a simple square pulse. By bringing $P$ and ${\mathcal{P}}$ closer via shaping the pulses, we can enable significantly better mode-parameter estimation. To this end, the elimination targets are:

(1) Cross Mode Coupling (CMC) error: Note an important difference between (2) and (3) is that former captures the coupling of internal ion states to modes $p\neq p^{*}$, an effect hereafter referred to as the CMC, whereas the latter does not. Specifically, the leading-order CMC error may be quantified as

$$\theta^{\rm CMC}_{p}=\Theta_{p}^{(1)}:=\int_{0}^{\tau}g(t)e^{i\omega_{p}t}dt\quad\quad(p\neq p^{*}),$$

(4)

where $\Theta_{p}^{(1)}$ is the first-order Magnus integral, since the first-order Magnus-approximated evolution operator $\hat{U}^{(1)}$ for (2) is

$$\hat{U}^{(1)}:=\exp(\hat{\Omega}_{1})=\exp\left(\sum_{p=0}^{N^{\prime}-1}\eta_{p}\Theta_{p}^{(1)}\hat{\sigma}^{+}\hat{a}_{p}^{\dagger}-h.c.\right)$$

(5)

and we used $\hat{\Omega}_{1}$ to denote the first-order Magnus term. Note this approximation is valid when $|\eta_{p^{*}}g(t)|$ is small compared to $1/\tau$. By eliminating all $\theta_{p}^{\rm CMC}$ defined in (4) using pulse shaping, the first-order CMC effect may be suppressed, leaving the second-order Magnus term

$$\hat{\Omega}_{2}=\sum_{p,p^{\prime}=0}^{N^{\prime}-1}\eta_{p}\eta_{p^{\prime}}\Theta_{p,p^{\prime}}^{(2)}\left(\hat{\sigma}_{z}\hat{a}^{\dagger}_{p}\hat{a}_{p^{\prime}}+\hat{\sigma}_{-}\hat{\sigma}_{+}\delta_{p,p^{\prime}}\right)-h.c.,$$

(6)

where

$$\Theta_{p,p^{\prime}}^{(2)}=\frac{1}{2}\int_{0}^{\tau}dt_{1}\int_{0}^{t_{1}}dt_{2}g(t_{1})g^{*}(t_{2})\>e^{i\omega_{p}t_{1}}\>e^{-i\omega_{p^{\prime}}t_{2}}$$

(7)

is the second-order Magnus integral, as the post-suppression leading-order CMC error.

(2) Detuning error: In our model Hamiltonian, $\omega_{p}$ are considered static. However, experimental imperfections, either from faulty characterization or drifts during the experiment Maksymov et al. (2022), give rise to inexactness or detuning of the mode frequencies $\omega_{p}$. This manifests in our model as

	$\displaystyle\theta^{\rm{det}}_{p}(\delta_{p})$	$\displaystyle:=\int_{0}^{\tau}e^{i(\omega_{p}+\delta_{p})t}g(t)dt-\int_{0}^{\tau}e^{i\omega_{p}t}g(t)dt$
		$\displaystyle=\sum_{\kappa=1}^{\infty}\frac{\delta_{p}^{\kappa}}{\kappa!}\frac{\partial^{\kappa}}{\partial\omega_{p}^{\kappa}}\Theta_{p}^{(1)},$		(8)

which is the difference incurred in the first-order Magnus integral due to the detuning $\omega_{p}\rightarrow\omega_{p}+\delta_{p}$. For the target mode $p^{*}$, the difference leads to off resonance of the BSB transition, which directly causes error in estimating $\eta_{j,p^{*}}$. For the non-target modes, the detuning leads to changes in the CMC terms in (4), hence imperfect removal of the CMC effects via pulse shaping, which ultimately results in the characterization error for $\eta_{j,p^{*}}$. Since the precise value of $\theta_{p}^{\rm{det}}$ cannot be learned when $\delta_{p}$ is not known precisely, we aim to suppress $\theta_{p}^{\rm{det}}(\delta_{p})$ via minimizing the first $K$ leading-order derivatives ($\kappa\leq K$) of $\Theta_{p}^{(1)}$ within the Taylor expansion terms of (8).

To recap, we emphasize that the quantities (4) and (8) are in general nonzero for the non-modulated pulse of the form $g(t)=\bar{A}e^{-i\mu t+i\phi}$. In our work, we deliberately design $g(t)$ to remove (4) and (8), such that the accuracy of mode characterization is significantly improved. Note that the residual errors in the characterization would then be due to the contributions from the second- and higher-order Magnus terms as mentioned earlier, in addition to other high-order Hamiltonian terms unaccounted for in our model. Indeed, an example may be the aforementioned DW effect. In the rest of our paper, we therefore neglect the DW effect and focus on the larger, leading-order effects such as the CMC and detuning error that we actively silence using our pulse shaping.

In this section, we present pulse-shaping techniques for suppressing the CMC and the detuning errors by removing (4) and (8), respectively. To summarize, we aim to obtain the following goals:

1. Goal 1

   Suppress the CMC by imposing $\Theta_{p}^{(1)}=0$ for all $p\neq p^{*}$;

2. Goal 2

   Maximize $|\Theta_{p^{*}}^{(1)}|$ while keeping the average Rabi frequency of the pulse the same, such that maximal change in the qubit population (observable in our experiment) is reached with a reasonable power requirement;

3. Goal 3

   Achieve Goal 1 and Goal 2 and estimate $\eta_{j,p^{*}}$ well even in the presence of uncertainties in $\omega_{p}$’s by stabilizing the conditions with respect to changes in $\omega_{p}$’s.

Briefly, to eliminate the CMC error, we first assume $\omega_{p}$’s are fully known, a condition which we relax later when introducing the detuning errors. Our Goal 1 is then to obtain $\theta_{p}^{\rm CMC}=0$ in (4) for all $p\neq p^{*}$. This way, the non-target modes are decoupled from the illuminated ion up to the first order.

In principle, a pulse that achieves Goal 1 can, for example, be far detuned from all of the mode frequencies. This would result in a minimal target-mode signal (qubit population) or, said differently, a pulse solution with a prohibitively large power. To address this, we maximize the target mode’s response (Goal 2) by maximizing $|\Theta_{p^{*}}^{(1)}|$ while keeping the pulse power requirement constant, to be more concretely defined later.

Relaxing now the perfect knowledge assumption for $\omega_{p}$’s, we notice that the precise integral values $\Theta_{p}^{(1)}$ used to achieve Goal 1 and Goal 2, for both $p{=}p^{*}$ and $p{\neq}p^{*}$, will change if a detuning error occurs. To prevent a significant change, we can stabilize all $\Theta_{p}^{(1)}$’s with respect to their respective mode frequencies $\omega_{p}$’s such that $\theta_{p}^{\rm det}$ in (8) is suppressed to (near) zero for a small range of detuning, with a varying degree of stabilization $K$. In other words, for Goal 3, we can minimize the derivatives of $\Theta_{p}^{(1)}$ with respect to $\omega_{p}$ for all modes $p$.

In what follows, we focus on achieving all three goals stated above by appropriately shaping $g(t)$. To start, we expand $g(t)$ using a Fourier basis, i.e.,

$$g(t)=\sum_{n}A^{\prime}_{n}e^{-i\frac{2\pi n}{\tau}t+i\phi_{n}}=\sum_{n}A_{n}e^{-i\frac{2\pi n}{\tau}t},$$

(9)

where $n$ indexes the Fourier basis and $A_{n}:=A_{n}^{\prime}e^{i\phi_{n}}$ is a complex coefficient of each Fourier component. While a complete Fourier basis comprising an infinite number of basis components may be considered, in practice, a Fourier basis consisting of $N_{\rm{basis}}$ components, concentrated around motional mode frequencies $\omega_{p}$’s, suffices. The average Rabi frequency $\bar{A}$ of pulse $g(t)$ is then defined as $\bar{A}:=\sqrt{\sum_{n}|A_{n}|^{2}}$. Inserting (9) into $\Theta_{p}^{(1)}$ defined in (4), we obtain

$\displaystyle\Theta^{(1)}_{p}=\sum_{n}A_{n}\int_{0}^{\tau}\exp\left\{i\left(\omega_{p}-\frac{2\pi n}{\tau}\right)t\right\}dt,$

(10)

which we can then summarize into a $N^{\prime}\times 1$ column vector $\vec{\Theta}^{(1)}$. Extracting the Rabi frequency coefficients $A_{n}$ from the above expression and summarizing again as a column vector $\vec{A}$, we can succinctly express $\vec{\Theta}^{(1)}$ as

$\displaystyle\vec{\Theta}^{(1)}=\mathbf{M}\vec{A},$

(11)

where $\mathbf{M}$ is a matrix of dimension $N^{\prime}\times N_{\rm{basis}}$ with the matrix elements

$\displaystyle M_{p,n}=\int_{0}^{\tau}\exp\left\{i\left(\omega_{p}-\frac{2\pi n}{\tau}\right)t\right\}dt.$

(12)

The goal here becomes then to determine the appropriate Fourier expansion coefficients $\vec{A}$ that achieve the three goals listed above.

Step 1 (Goal 1): CMC suppression – Here, we aim to ensure the inner product between each row of $\mathbf{M}$ and $\vec{A}$ is zero for all $p{\neq}p^{*}$, such that $\theta_{p}^{\text{CMC}}{=}0,\,\,\forall p{\neq}p^{*}$. To achieve this, we remove the $p^{*}$-th row vector from $\mathbf{M}$, denote this matrix as $\mathbf{M}^{\prime}$, and find the null space of $\mathbf{M}^{\prime}$. We characterize this space $\boldsymbol{\mathcal{S}}^{(1)}$ by a spanning set of vectors $\mathcal{S}^{(1)}$, and denote $\mathbf{N}$ as the matrix whose columns consist of vectors in $\mathcal{S}^{(1)}$. These nullspace vectors then satisfy the CMC nulling conditions of all $p\neq p^{*}$ and serve as the spanning vectors of the space over which we maximize the signal of the target mode $p=p^{*}$. The size of the nullspace is the number of basis functions, subtracted the number of nulling conditions applied $N^{\prime}-1$. The null-space dimension marks the degrees of freedom we have left to work with for the signal maximization, of which there needs to be a sufficient number to warrant convergence in the pulse shape. This implies that the pulse length has to be sufficiently long such that the frequency resolution $2\pi n/\tau$ is small enough to generate sufficiently many frequency basis components.

Step 2 (Goal 2): Maximizing signal strength – To maximize the signal strength of the target mode $p^{*}$, we first minimize the average power of the pulse to achieve $|\Theta_{p^{*}}^{(1)}|=1$. We then scale the obtained, minimal-power pulse solution vector by Rabi-frequency scaling factor $\alpha$, such that $|\Theta_{p^{*}}^{(1)}|$ reaches a desired value, in this case $\alpha$. More specifically, our starting point is the projection of the $p^{*}$-th row vector $\vec{v}^{\dagger}$ of $\mathbf{M}$ onto the nullspace $\mathbf{N}$. We then find the eigenvector $\vec{\xi}$ with the largest eigenvalue $\lambda_{\rm max}>0$ of the positive-definite matrix $\mathbf{N}^{\dagger}\vec{v}\vec{v}^{\dagger}\mathbf{N}$. The normalized pulse solution is obtained by projecting $\vec{\xi}$ onto the nullspace $\mathbf{N}$ then normalizing by $\sqrt{\lambda_{\rm max}}$. After rescaling by a factor of $\alpha$, we obtain the pulse solution $\vec{A}$, given by

$\displaystyle\vec{A}=\frac{\alpha}{\sqrt{\lambda_{\rm max}}}\mathbf{N}\vec{\xi}.$

(13)

It is worth mentioning that our choice of $\alpha$ is restricted by the conditions of the RWA approximations. As a result, we mainly operate in the small $\alpha$ regime.

(Optional) Modified Step 1 (Goal 3): Stabilization against detuning – We now admit that we have inaccuracies in $\omega_{p}$’s such that $\omega_{p}\rightarrow\omega_{p}+\delta_{p}$, $\forall p$. The resulting error contribution $\theta_{p}^{\rm det}$, given in (8), can be suppressed by adding the following constraints to the null-space conditions:

$$\frac{\partial^{\kappa}}{\partial\omega_{p}^{\kappa}}\Theta_{p}^{(1)}=0\quad\quad\forall p,\>\kappa\in\{1,..,K\},$$

(14)

where $K$ is the moment of stabilization, which corresponds to the highest order of derivatives that are nulled.

Specifically, these constraints are appended as additional $KN^{\prime}$ rows to the matrix $\mathbf{M}^{\prime}$, where each row contains the nulling condition for each pair of $p$ and $\kappa$. As $\mathbf{M}^{\prime}$ originally has $N^{\prime}-1$ rows, the appended matrix entries in the expanded $\mathbf{M}^{\prime}$ are given by

$\displaystyle M^{\prime}_{\kappa N^{\prime}-1+p,n}=\frac{\partial^{\kappa}}{\partial\omega_{p}^{\kappa}}\int_{0}^{\tau}\exp\left\{i(\omega_{p}-2\pi n/\tau)t\right\}dt,$

(15)

where $p=1,..,N^{\prime}$ and $\kappa=1,..,K$. Then, we find the set of null space vectors $\mathcal{S}^{(1)}$ corresponding to the expanded $\mathbf{M}^{\prime}$ and proceed to perform the same signal maximization procedure as described in Step 2.

We note in passing that performing stabilization against detuning as described in Modified Step 1 is optional. Should the rough $\omega_{p}$ values known a priori to the characterization be sufficiently accurate for the purpose of characterization, the need for these constraints are obviated.

Remarks – Due to its role in the rescaling of $|\Theta_{p^{*}}^{(1)}|$, hereafter, we use $\alpha$ interchangeably as the absolute value of the first-order Magnus integral for mode $p^{*}$ for both shaped and square pulses. The rescaled pulse roughly induces a qubit population inversion $\mathcal{P}\approx\sin^{2}(\eta_{p^{*}}|\Theta_{p^{*}}^{(1)}|)$ (see Appendix B.3 for details) and, in the perturbative regime where $|\eta_{p^{*}}\Theta_{p^{*}}^{(1)}|\ll 1$, $\mathcal{P}\approx|\eta_{p^{*}}\Theta_{p^{*}}^{(1)}|^{2}=(\eta_{p^{*}}\alpha)^{2}$. The pulse, if it happens to be that its Fourier coefficients for a tight band of frequencies near $\omega_{p^{*}}$ dominate in their modulus, would roughly have an average power of $\bar{A}\approx\alpha/\tau$ [see (10)].

In this section, we demonstrate the viability of our pulse-shaping techniques proposed via numerical simulations of a three-ion chain. To benchmark the performance of our pulse-shaping characterization tools as compared to their traditional square-pulse counterparts, we simulate the bright-state populations $P$ and $\mathcal{P}$ induced by the multi-mode Hamiltonian (2) and the single-mode Hamiltonian (3), respectively, – see Appendix A of Ref. Kang et al. (2023a) for the simulation implementation detail – where the initial qubit state is $\ket{0}$ and all motional modes are initially in the ground state. We then use the fractional difference between the induced qubit populations by the two models $\mathcal{E}:=|P-\mathcal{P}|/\mathcal{P}$, hereafter referred frequently as the fractional population error, resulting from the CMC and detuning error, as a proxy for the characterization error.

It is worth noting that since $P\approx\sin^{2}(\eta_{p^{*}}\alpha)$ (see Appendix B.3), $P\approx(\eta_{p^{*}}\alpha)^{2}$ in the small-population regime, so the fractional error in $P$ in relation to the fractional error in $\eta_{p^{*}}$ is given by $\frac{dP}{P}\approx\frac{2d\eta_{p^{*}}}{\eta_{p^{*}}}$. Thus, in the presence of fitting errors or experimental imperfections, $\mathcal{E}$ serves as a lower bound to the fractional characterization uncertainty of $\eta_{p^{*}}$. Computing the proxy error $\mathcal{E}$ for both shaped and square pulses can hence help us reveal the parameter regime in which shaped pulses offer an advantage over square pulses.

Specifically, to see how the parameter regime may be found, we can consider the leading-order errors that remain after pulse shaping and contribute significantly to ${\mathcal{E}}$, which can be broken down as follows. First, residual CMC errors of order $O(\eta^{2})$ or higher remain due to the second- or higher-order Magnus terms. Second, detuning errors arise from higher-order contributions of the first-order Magnus integral $\frac{\delta_{p}^{\kappa}}{\kappa!}\frac{\partial^{\kappa}}{\partial\omega_{p}^{\kappa}}\Theta_{p}^{(1)}$ ($\kappa>K$), where $K$ is the moment of stabilization. Third, additional detuning errors arise from non-target higher-order Magnus integrals. Thus, the dependence of $\mathcal{E}$ to various parameters such as pulse power, pulse duration, and detuning can be predicted based on the residual error sources listed above. Indeed, we aim to identify parameter choices that result in smaller $\mathcal{E}$ for shaped pulses compared to that for square pulses by carefully examining each dependence.

We note in passing that, for a linear chain of three ions indexed by $j\in\{0,1,2\}$ (labeled from left to right) with three modes indexed by $p\in\{0,1,2\}$ (labeled in the ascending order of mode frequencies), we focus throughout this section on ion $j=2$, probing the largest-frequency mode $p^{*}=2$, as a concrete example; This choice is arbitrary and other ion-mode pairs can be considered straightforwardly in our procedures. Further, we use a realistic set of values of $\omega_{p}$ and $\eta_{j,p}$, reported in Appendix A, for our simulated experiments (three-mode Hamiltonian) throughout this section.

Figure 1 shows an example square pulse and shaped pulses with various moments of stabilization $K$, for the aforementioned example choice of the ion and the mode. The pulse components in the frequency domain are mostly concentrated near the target-mode frequency $\omega_{p^{*}}{=}\omega_{2}$. The small wiggles around the non-target mode frequencies $\omega_{0}$ and $\omega_{1}$ can be interpreted as signatures of active cancellation of the CMC error. Similar features are observed for other choices of ion-mode pairs (not shown).

The rest of the section is organized as follows. We first study the $\mathcal{E}$ scaling of square and shaped pulses without stabilization due to the CMC effect in Sec. V.1. We then study the $\mathcal{E}$ scaling of square and shaped pulses with various degrees of stabilization in the presence of both errors, due to CMC and detuning, in Sec. V.2. Finally, we compare $\mathcal{E}$ between various shaped and square pulses across all parameter regimes in Sec. V.3.

We devised a pulse-shaping scheme for high-accuracy characterization of the motional-mode parameters of trapped-ion quantum computers. We have extensively tested the effectiveness of our scheme in minimizing the CMC and detuning errors via simulations on a three-ion chain. Our results showed an improvement in the accuracy of $\eta_{j,p}$ estimation, even when the mode frequencies were not well known.

Our pulse-shaping scheme helps address one of the many engineering challenges for building large trapped-ion quantum computers. As the number of qubits within a trap increases, an accurate and efficient mode characterization in the presence of other motional modes becomes increasingly challenging. By leveraging the degrees of freedom in pulse design for these experiments, our methodology offers one potential path forward to overcoming this bottleneck, namely, accurately probing mode parameters in the presence of model violations and experimental imperfections. One of the exciting future directions may be to explore parallel pulse shaping, where multiple ions are addressed simultaneously to characterize mode parameters more efficiently.

The accurate and efficient characterization of these mode parameters holds the key to high-fidelity trapped-ion operations, such as parallel entangling gates Bentley et al. (2020); Grzesiak et al. (2020). Briefly, these gates are implemented by illuminating each ion by custom-designed pulse shapes such that the induced interactions between the internal states of the ion qubits match those intended as quantum gates acting over multiple pairs of qubits simultaneously. Such interactions are, just as in a single two-qubit gate between two ion qubits, mediated by motional modes of the ion chain. As such, the parallel entangling gates intricately depends on the precise knowledge of the mode parameters for their high-fidelity implementation. Provided that these parallel gates have been proven to be versatile in enabling efficient quantum computing over a wide class of quantum programs Grzesiak et al. (2022); Bravyi et al. (2022), the improvement achieved in characterizing mode parameters enabled by our methodology is expected to play a crucial role in realizing efficient quantum computation.

Our methodology of applying pulse shaping and other optimal control techniques for characterization tasks may readily be adopted to be applicable for different quantum computer architectures as well. Systems that rely on spectroscopic system parameter characterization like trapped ions can likely benefit from pulse shaping in their characterization by effectively eliminating significant model violations. We hope our work, which opens a large trade space of system-level optimization for exploration and exploitation, will help boost the effort of building a large-scale quantum computer.

Q.L. and M.K. thank Kenneth R. Brown for helpful suggestions. M.K. was supported by the Army Research Office, W911NF-21-1-0005.