Sculpting bosonic states with arithmetic subtractions

Continuous-variable (CV) encoding allows information to be processed compactly and efficiently on quantum processors. Recently developed techniques such as controlled beam-splitter operations and the near deterministic phonon subtractions make trapped ion systems attractive for exploring CV quantum computing. Here we propose a probabilistic scheme based on the boson sculpting technique for generating multipartite highly entangled states of motional modes of trapped ion systems. We also investigate the effects of decoherence on the fidelity of the generated state by performing numerical simulations with realistic noise parameters. Our work is a step towards generating multipartite continuous-variable entanglement.


I. INTRODUCTION
Quantum entanglement is a property of a compound system that possesses non-classical correlation between its subsystems. The ability to prepare highly entangled states is crucial for quantum computation. For physical platforms that employ discrete variables, the entanglement between qubits are typically created using shortrange interactions induced by bosonic modes. In particular, with the collective motional modes of trapped ions as the quantum bus, the internal degrees of freedom of ions have been entangled with a fidelity significantly above the threshold required for fault-tolerant quantum computation [1][2][3]. However, with a larger number of qubits, implementing the full control necessary for entangling operations remains challenging due to several technical problems, such as crosstalk, heating, and the overhead of addressing individual qubits.
Alternatively, classical information can be encoded into the eigenstates of continuous-valued operators, such as the motion of trapped ions, the quadratures of an electromagnetic mode, and the spin variables of an atomic ensemble [4][5][6]. In such bosonic systems, a large dimension of the Hilbert space is typically available for the encoding. As less physical resources are required, this makes the quantum computation more efficient.
The bosonic system we study in this paper is the trapped-ion system. The motion of a linearly trapped ion chain is a well-controlled bosonic mode, which allows for deterministic preparation of single-phonon states [7]. There have also been several experiments that report nonlinear gates operating on trapped-ion phonons with good fidelity [8,9], and high-fidelity state reconstruction of the motional state is also possible [10].
Despite their advantages, it can be more complicated to generate entanglement with bosonic modes, though there have been some early experimental attempts to generate bipartite entanglement [11][12][13]. Here, we present a * These authors contributed equally to this work. scheme based on arithmetic phonon subtraction to prepare a cat state on the motional modes of 4 ions. The proposed method can be generalized to a Hilbert space a higher dimension and more motional modes.
Cat states are evenly populated superpositions of maximally distinguishable states [14], also referred to as Greenberger-Horne-Zeilinger (GHZ) states. Because they contain genuine multipartite correlations, they are of particular interest for future quantum technologies [7,15]. This type of entanglement is considered a universal resource for quantum computing [16][17][18][19], quantum communication [20,21], and also for testing the foundations of quantum physics [22]. Cat states are challenging to prepare because of their sensitivity to decoherence. They induce the so-called super-decoherence which can be used as benchmark for robust quantum control [23][24][25][26][27].
In this paper, we propose an experiment for entangling the collective motional modes of 4 trapped ions. The same protocol can be extended to systems with a higher number of ions. The main idea of our approach is to prepare a single phonon in each motional mode, rotate the basis, and subtract half of the phonons. This results in a final state that carries genuine multipartite correlations between the motional modes. Because particle subtraction is the key step for revealing these mode correlations, we refer to this approach as a sculpting scheme, which was originally proposed for multimode photonic platforms [28]. Since it is not possible to remove particles from a vacuum, sculpting schemes are probabilistic, with the success probability dependent on the vacuum component of the state prior to the subtraction process.
As quantum entanglement is regarded as a vital resource for quantum technologies, finding new ways for its extraction is of natural interest. Moreover, for manybody systems composed of indistinguishable particles, the intrinsic correlations due the symmetrization constraints are the subject of a rapidly growing interest for their potential applications as a quantum resource [29][30][31][32]. In the sculpting scheme, these intrinsic correlations are consumed in order to create entanglement between the bosonic modes. Here, we show that the sculpting scheme can also be implemented with trapped ions.
This article is organized as follows: in Section II, we start by presenting the basic idea of using subtraction of indistinguishable bosons for the creation of entanglement between bosonic modes. Then, in Section III, we discuss the necessary operations for adapting this scheme to trapped-ion platforms. In Section IV, we show how these operations can be used to generate entanglement from different initial states. In Section V, we discuss the numerical simulations of an experimental implementation assuming realistic conditions. In the last section, we give an outlook of this work and potential applications.

II. SCULPTING BOSONIC GHZ STATES
A. Creation of mode entanglement by subtraction Let us assume that we have 2n bosonic modes. We start by creating a single boson at each mode, corresponding to the symmetric state such thatâ † q is the creation operator acting on the q th mode and |∅ = |0 1 , · · · , 0 2n corresponds to the vacuum for all modes. Here we adopt the shorthand notation |n 1 , 0 2 , · · · , 0 2n which denotes having n particles at the first mode while all remaining modes are empty.
The initial symmetric state |sym 2n can be transformed into an entangled state by implementing n successive subtraction operations defined aŝ where each subtraction operation has the form After normalization, we obtain an entangled state of the form [28] Since our system consists of identical particles distributed over multiple modes, two types of quantum correlations are present: particle entanglement and mode entanglement [29][30][31][32]. The former is due to the exchange symmetry of the indistinguishable bosons, while the latter corresponds to correlations between different modes. In first quantization, the initial state reads such that σ is the set of all possible permutations of having one particle at each mode. Clearly, due to the exchange symmetry, the state above is highly correlated.
Since the entangled parties are indistinguishable, the correlations in |sym 2n are inaccessible. Through the transformation |sym 2n → |GHZ 2n , one can notice the following: (i) The amount of particle entanglement is reduced. This is because the number of particles decreases from 2n to n. Consequently, in first quantization, the set of possible permutations contracts as well. (ii) Mode entanglement of the GHZ -type is created. The choice of |sym 2n is important, as it has been shown that the initial state must have nonzero particle entanglement in order to extract accessible mode entanglement using only subtraction operations [30]. Some subtlety is involved in quantifying exactly what is exchanged between the two types of entanglement, which is briefly discussed in appendix A. In short, by performing n subtractions, some of the inaccessible particle entanglement is consumed in order to create accessible entanglement between bosonic modes [28].
In this paper, we shall rather focus on the so-called arithmetic operations [8] The arithmetic operations are referred to as "near deterministic" [8] because the only deviation from unitarity is due toŜ|0 = 0: that is, whileŜŜ † = I holds, one findŝ S †Ŝ = I − |0 0|. In particular, the subtraction S preserves the scalar product of all pairs of states (hence, the norm of all states) that do not have a vacuum component. In ion traps, arithmetic operations can be implemented without ancillary ions via adiabatic schemes [34][35][36]. If the adiabatic passage is performed slowly enough, and for a long enough time period, one can be certain that the populations have been completely transferred from each state |n to |n + 1 while maintaining the coherences [8].
In principle, the inclusion of arithmetic subtraction to the set of Gaussian operations performed in a trapped-ion system allows for universal state preparation [37]. However, there is no general method available to find the sequence of arithmetic subtraction and Gaussian operations required to prepare an arbitrary state. As such, we focus on reporting the exact state preparation of the |GHZ 4 state, which has a known use in quantum computation [38].

C. Bosonic sculpting with trapped ions
In trapped-ion systems, we can address the motional modes of the ions in two bases: the local basis, and the collective basis. The local basis, whose states and operators we will denote with a subscript l, corresponds to the motion of the l th ion. The collective basis, whose states and operators we will denote with a subscript c, refers to the normal modes of the collective motion of the ion chain.
In this proposal, we will show how to create the target state |GHZ 4 c , entangled in the collective basis, starting from an initial state prepared either in the local basis |sym 4 l or the collective one |sym 4 c . The gate sequence for both scenarios is shown in Fig. 1. The first scenario (subsection II C 1) requires fewer gates, but is harder to implement faithfully because of the phonon-hopping between the local modes. The second scenario (subsection II C 2) avoids this problem by working solely in the collective basis, but requires an implementation time that is almost thrice as long, and is thus more affected by noise.
The initial state is prepared in the local basis |sym 4 l , and the first gate to act on it is a beam-splitting operation between the 2 nd and 4 th collective modesB 2,4 ≡ B 2,4 (2λ − π 2 , − π 2 ). For the sake of convenience, we define here the following ion-trap beam-splitting convention: The implementation details of these beam-splitting gates will be discussed later (see subsection III C 3). Using the above definitions in (9) and (8), we can writẽ By expanding the product above, one can easily see that after subtracting one particle each from the 3 rd and 4 th collective modes, only two terms will survivê such that Comparing this to the target state, we obtain the ideal fidelity | ψ f |GHZ 4 c | 2 = 3/ √ 10 ≈ 0.949. With arithmetic subtractions, we do not obtain |GHZ 4 c with unit fidelity:Ŝ c,k does not give rise to the usual factors of √ n like the ordinary annihilation operatorâ c,k , so the two constituent states making up the superposition in |ψ f are not evenly populated. While this state is no longer separable, it is not yet maximally entangled. In order to increase the correlations between the motional modes of the ions, we need to evolve the state as Here, |g /|e refers to the internal spin state of the ion addressed by the Raman lasers-see subsection III A for more details. The transformation above can be implemented using the red sideband transition on the 2 nd collective mode, which will be discussed in subsection III B 2. As a final step, we post-select the ground state |g to get, after normalization, (14) All beam-splitting operations considered in this article are between collective modes. However, if the initial state is prepared in the local basis as |sym 4 l = a l,1âl,2âl,3âl,4 |∅ , one will need to use the relations (8) to appropriately compute the output state. If the initial state is prepared in the collective basis |sym 4 c = a c,1âc,2âc,3âc,4 |∅ , then the transformation in (9) can be used in a straightforward manner.   (9), the beam-splitting gateB4,2 has the parameters φ = −π/2 and θ = 2λ − π/2. (b) The initial state is prepared in the collective basis. All the gates Bp,q are 50-50 beam-splitting operations, i.e. φ = −π/2 and θ = π/2. The gateŜ c,k correspond to an arithmetic subtraction from the k th collective mode. The implementation of the arithmetic subtraction process involves post-selecting the ground state of the spin degree of freedom, as detailed in subsection III C 2. The remaining gates correspond to a red sideband transition, followed by another post-selection of the ground state, to obtain exactly the target state |φ4 c.

Scenario without individual addressing
In the second scenario, the ions are not addressed individually. The initial state is prepared in the collective basis |sym 4 c . First, let us define a sequence of 50-50 beam-splitting operations as follows Using the definition (9), one can easily verify that the operationB j,k ( π 2 ,− π 2 ) is a 50-50 beam-splitting gate between the j th and k th modes. If we apply sequence above to the initial state, we get This is exactly the state (10) that we previously obtained after the beam-splitting gateB 2,4 . Therefore, to reach the target state |GHZ 4 c , one has to perform the same sequence of operations as before (subtract from the 3 rd and 4 th modes, drive the red sideband transition of the 2 nd mode, and finally post-select the ground state of the internal degree of freedom). In general, this scheme can be extended to a |GHZ 2n state for any n, with and without the red sideband correction. We discuss the general case in appendix B. For this experimental proposal, we focus only on the simplest nontrivial case of n = 2.

Success probabilities
Due to the post-selective measurements, required for both the arithmetic subtraction and the corrective red sideband operation, this scheme is probabilistic. The success probabilities are 5/16 = 31.25% without the red sideband correction, and 1/8 = 12.5% with the red sideband correction. These probabilities are the same for the scenario with and without individual addressing.
In trapped-ion systems, the initial state |sym 2n is deterministically prepared from the motional ground state with the red sideband transition, as will be covered in section III B 2. The measurement of the internal state of the ion is also highly efficient, with a state detection fidelity of 0.999 [41]. Therefore, the success probability of preparing |GHZ 2n is determined solely by the success probability of the subtraction sequence and the red sideband correction, as reported above.

III. BOSONIC MODES OF TRAPPED IONS
A. Hamiltonian of the system For the concrete proposal, we refer to the experimental setup reported in [6,42,43]: a system of four ions in a linear Paul trap, with trap frequencies ω z ω x , ω y . We address the motional degree of freedom of the ions along the x direction and the internal spin state of the last (fourth) ion. The last ion is chosen since it can participate in all motional modes. The Hamiltonian of the system takes the form [40,44] Here, ν j andâ c,j are the frequency and annihilation operator of the collective mode j, ω 0 is the carrier transition frequency, andσ z = |e e| − |g g|, where |e and |g are the excited and ground states of the spin.
The system is driven by a pair of lasers that couple the internal states of the fourth ion to the collective motional modes via Raman transitions. We denote the laser frequencies as ω L,1 and ω L,2 , and their phases as φ 1 and φ 2 . In the rotating frame of the free HamiltonianĤ 0 , and after taking the rotating-wave approximation, the ion-laser coupling is governed by the interaction Hamiltonian [45] wherê Here, {â c,j } 4 j=1 are annihilation operators acting on the collective motional modes, while {η j } 4 j=1 are the Lamb-Dicke parameters.
The terms in each brace in equations (19)(20)(21)(22)(23) correspond to a resonant frequency. They be addressed by adjusting the effective laser detunings δ 1 := ω L,1 − ω 0 and δ 2 := ω L,2 − ω 0 such that the corresponding terms in the Hamiltonian become time independent, while the contributions of the rapidly oscillating off-resonant terms become negligible.
By addressing different resonant frequencies, the various quantum operations required for the boson sculpting scheme can be performed. We list these quantum operations, and the laser parameters required to perform them, in Table I.

B. Basic operations
In the following sections, the coupling strengths of the Raman lasers are taken to be equal and time dependent: that is, g 1 = g 2 = g(t). Furthermore, all relevant Hamiltonians will be found to be in the formĤ(t) = g(t)Ĝ for some time-independent operatorĜ. The time evolution operator of a system evolving under such a Hamiltonian is given byÛ (t) = exp(−i t 0 dt g(t )Ĝ).

Carrier Transition
The carrier transition is zeroth order with respect to the Lamb Dicke parameter, which can be performed by setting the laser detunings to δ 1 = δ 2 = 0, and phases to The time evolution operator of this Hamiltonian, U carr (θ, φ), with θ = 2 t 0 dt g(t ), is immediately recognizable as a rotation of the spin statê This transition is required for the preparation of motional states in conjunction with the red sideband transitions, and for setting the spin to the correct state to perform operations on the motional degree of freedom.

Red Sideband Transition
The red sideband transition is a first order operation that can be performed by setting the laser detunings to δ 1 = δ 2 = −ν j , and phases to φ The time evolution operator U rsb,j (θ, φ), with θ = 2η j t 0 dt g(t ), results in Rabi oscillations between the states |e |n c,j ↔ |g |n + 1 c,j , where At this point, one can easily verify the transformation (13) by substituting the values θ = 2π/3 and φ = π/2 into (27)(28).
Also, this transition can be used to prepare the initial state required for the boson sculpting scheme. Consider the state |g |0 c,j . If a carrier transition is performed, followed by a red sideband transition on mode j, with θ = π, φ = π/2 for both transitions, the state undergoes the evolution Repeating this sequence of transitions for all modes j, we would have

Blue Sideband Transition
The blue sideband transition is also a first order operation, which can be performed by setting δ 1 = δ 2 = ν j and φ 1 = φ 2 = −φ − π 2 , such that The time evolution U bsb,j (θ, φ), with θ = 2η j t 0 dt g(t ), carries out the transformation U bsb,j (θ, φ)|e |n + 1 c,j = cos( θ √ )|g |n c,j , Generally, operations involving the red sideband transition can also be performed with the blue sideband by first performing a π-rotation on the spin. In this paper, we address the red sideband instead of the blue sideband wherever required.

C. Composite operations
In the previous subsection, we presented the 3 basic transitions necessary to manipulate the internal and motional degrees of freedom of the trapped ions. In the following, we will show how these transitions can be used to implement Gaussian gates such as displacement and beam-splitting operations. While the former will be necessary for the tomography of the system, the latter will be useful for preparing the ions in the appropriate basis. Also, we can implement non-Gaussian gates, namely, arithmetic operations, which will be used for subtracting phonons from our system.

Displacement Operation
This is a first order operation performed by driving both sideband transitions simultaneously with By setting the spin state of the system to |− := 1 √ 2 (|e − |g ) with the carrier transition, the time evolution operator corresponds to a displacement opera-tionD j (θ, φ) := exp θ(â † c,j e iφ −â c,j e −iφ ) , with θ = ηj 2 t 0 dt g(t ), which carries out the transformation While the boson sculpting scheme does not require the use of the displacement operator, it can be useful for tomography when used with the parity gate.

Arithmetic subtractions of phonons
The arithmetic subtraction operator on mode k can be performed by adiabatically driving the red sideband transition given in (26). This is done by slowly varying the laser detuning and coupling strengths over a time interval t ∈ [0, τ ] as [8] δ 1 (t) = δ 2 (t) = −ν k + ∆ 0 cos( πt τ ) , where the detuning ∆ 0 = 1 2 √ n max + 1η k g 0 depends on n max , the maximum number of phonons to be subtracted. When τ is large, the adiabatic theorem states that the eigenstates of the initial Hamiltonian evolves to the corresponding eigenstates of the final Hamiltonian with the same eigenvalue [46]. In this case, this results in an adiabatic transfer of the states from |g |n c,k → |e |n − 1 c,k , given that the adiabatic condition τ 1/gη k is met [8]. A final carrier transition is performed to reset the spin to |g |n − 1 c,k . As the red sideband transition does not affect the state |g |0 c,k , this final transition excites the spin to the state |e |0 c,k . A projective measurement on the spin thus removes the vacuum component of the state, completing the arithmetic subtraction process.

Beam-Splitting operations between the collective modes
In recent years, significant progress have been made studying the motional degree of freedom of trapped ions [47]. This has led to the demonstration of new tools and techniques that allow for better control over such systems. For instance, due to Coulomb interactions between ions, a beam-splitter-like coupling between the motional modes of two trapped ions has been experimentally demonstrated [33]. This has led to many applications, such as the implementation of the controlled-SWAP gate for machine learning algorithms [6,43] and the study of quantum walks using phonons [48].
The beam-splitter transformation is a second order operation that can be performed by setting δ 1 = −δ 2 = ν j − ν k and φ 1 = −φ 2 = π − φ, so that H B,j,k (t) = σ x g(t)η j η k 2 â † c,jâ c,k e iφ +â c,jâ † c,k e −iφ . (33) By setting the state of the spin to |− = 1 √ 2 (|e − |g ), the time evolution is exactly the unitarŷ with θ = η j η k t 0 dt g(t ). This carries out the beamsplitting transformationŝ As we will see in the following section, the beam-splitting operation is vital for preparing the initial state in the boson sculpting scheme. Furthermore, the parity gate can be performed by choosing different values for the φ 1 and φ 2 , as detailed in [42]. Together with the displacement operation covered in the previous section, a direct measurement of the joint Wigner function can be performed [49].

IV. STUDY OF THE EFFECTS OF HEATING AND DECOHERENCE
Now that we have defined the sequence of operations for our scheme, we simulate its experimental implementation while considering realistic noise conditions in this section. Here, we assume two sources of external per-turbation: the damping of the motional modes due to fluctuations of the trap frequency, and the heating of the motional modes. We model the first effect by an interaction with a phase damping reservoir, and the second with an amplitude damping reservoir. Then, the time evolution of the quantum state is computed via the master equation [50] where {·, ·} is the anticommutator, γ r and κ r are the decay rates due to the coupling to the amplitude-and phase-damping reservoirs of the r th motional mode respectively, andn = 10 6 is the average number of phonons in the reservoir. For a given on-resonant transition, the angle of rotation is proportional to the pulse area τ 0 dtg(t). Typically, to mitigate undesired dynamics involving other transitions, the applied pulse shapes have smooth rising and falling edges. In our numerical simulations, we chose a soft-edged square pulse of the form The area of the pulse is τ 0 dtg(t; g 0 , τ, t r ) = g 0 (τ − t r ), where g 0 is the amplitude and t r the rise time of the pulse. The pulse amplitude depends on the power output of the laser, where we used g 0 = π/0.004 [43]. Meanwhile, a longer rise time reduces the noise due to off-resonant transitions at the cost of a longer pulse time. Unless otherwise stated, we use t r /τ = 0.125 in our numerical simulations.

A. Numerical simulation of the gates
Before considering the entire sequence of operations given in Fig. 1, which gives rise to the target state,Û  we shall study the effects of noise on each gate separately. In these numerical simulations, contributions of the off-resonant terms are included by considering the full Hamiltonian (18). Heating and decoherence effects are included using Eq. (35), with decay rates nγ 01 = 15 phonons/s,nγ 0 (2,3,4) = 0.675 phonons/s, and κ 0 = 0.075 ms −1 . To implement each gate, the laser detunings are set according to Table I. In Fig. 2, we plot the infidelity of the sideband transition while varying the pulse length. We address the second collective motional modes in order to get the transformation (28), i.e. |g |1 2 c → |e |0 2 c . For a duration up to a π-pulse, this gate exhibits a robust behaviour with infidelities bellow 10 −3 . The blue dot in Fig. 2 corresponds to the transformation (13). In the inset, we plot the the population transfer from |g |1 2 c to |e |0 2 c .
In Fig. 3, we plot the infidelities of the beam-splitting gate between several pairs of modes. Note that the pairs chosen here are the same pairs discussed in subsection II C. Addressing these pairs is needed for the implementation of the entangling scheme. Namely, the diamond dot corresponds to the single beam-splitting transformation required for implementing the scenario with individual addressing II C 1. On the other hand, the circle dots correspond to the 4 beam-splitting transformations required for the implementation of the scenario without individual addressing II C 2. Since the angle of rotation of any beam-splitting transformation is proportional to the pulse area, larger angles will require longer gate times and hence more errors will accumulate. Consequently, one might predict that the scenario with individual addressing will be more robust against noise. such a prediction will be confirmed in the following subsection. In the inset, we illustrate an example of a beam-splitting transformation between the 1 st and 2 nd modes with a single phonon initially in the 1 st mode. After a π-pulse, the phonon is transferred to the 2 nd with a ∼ 92% fidelity.
In Fig. 4, we plot the infidelity of the arithmetic subtraction gate as a function of duration of the adiabatic passage. For this gate, we can see that there is a an optimal time duration for which the infidelity is minimized. If the pulse length is shorter, the fidelity of the gate suffers because the adiabatic passage to too fast. On the other hand, if the duration is longer, the effects of noise become dominant.
B. Numerical simulation of the sculpting scheme In Figure 5, we can see that the fidelities in the scenario with individual addressing are considerably more robust against the noise. This can be explained by the fact that the run time for the scenario without individual addressing is almost 3 times longer that the other scenario (cf. Table II). For the scenario with individual addressing, the computed fidelity is 0.997 for the case without noise (bottom left of the color map). For the extreme case with only decoherence, κ = 0 and γ = 2γ 0 , it is found to be 0.917. For the case with only heating, κ = 2κ 0 and γ = 0, the fidelity is equal to 0.846. In the upper right corner, where the system is coupled to both baths with κ = 2κ 0 and γ = 2γ 0 , the fidelity reads 0.778. In the scenario without individual addressing, following the same order, the fidelities at the corners of the color map are 0.983, 0.733, 0.445, and 0.337.
In the case where the subtraction operation is performed with the adiabatic blue sideband, we instead obtain the fidelities 0.923, 0.859, 0.811, and 0.754 at the corners of the color map with individual addressing; and 0.921, 0.693, 0.434, and 0.331 without individual addressing. When using the blue sideband instead of the red, the behaviour of the system in the presence of noise is similar with slightly lower fidelities.
In subfigures (c) and (d), we plot the entries of the density matrix of the final state as a function of the coupling strength to the reservoirs. Overall, the subtraction scheme is more sensitive to the damping in both scenarios. However, the difference in sensitivity is relatively lower for the scenario with individual addressing.
In the cases without noise, the fidelities of the final state are not perfect despite the unitary dynamics of our system. This is mainly due to the contributions of the off-resonant terms of the Hamiltonian, as discussed in subsection III A. By comparing the isolated fidelities to the isolated run times in Table II, one might expect that the longer the run time of a gate, the more it accumulates such contributions (errors), and the less its isolated fidelity. While this might be intuitive, the cases of B 3,4 and B 2,4 are an example that such an intuition is not always correct. In this example, B 2,4 runs for a shorter time, but the isolated fidelity of B 3,4 appears to be higher. This can be explained by the fact that each gate corresponds to a different frequency detuning, ergo the con- c,4 |∅ . The cases where the system is coupled to both reservoirs (κ/κ0 = γ/γ0 = ξ) is labelled by the solid lines/circles, to just the heating reservoir (γ/γ0 = ξ, κ = 0) by the dashed lines/squares, and to just the damping reservoir (κ/κ0 = ξ, γ = 0) by the dotted lines/triangles. tributions of the remaining terms of the Hamiltonian are always different. These contributions can be suppressed by optimizing the shape of the laser pulses. While we used a soft-edged square pulse for all the operations in this work, identifying the most optimal pulse shape remain an open problem. This will allow for achieving even higher fidelities of the final state. From Table (II), one might notice that the isolated fidelities of the subtraction gates are not the same. The first subtraction from mode 3 has an isolated fidelity lower than that of the second subtraction from mode 4. This is because the gate times were optimized for the highest accumulated fidelity after both subtractions, as we are only interested in the resulting state given in equation (11). Therefore, this leads to a lower isolated fidelity of the intermediate state between the two subtractions.
The tomography of such states can be performed in two ways. The first method is to reconstruct the density matrix of the final state by performing projective measurements (see the supplementary material of [42]). Using the red sideband transition, the state encoded in the motional degree of freedom can be transferred to the internal degree of freedom, and beam-splitting operations such as B 12 (θ, φ) and B 34 (θ, φ) can be used to perform the necessary rotations. The second method is to reconstruct the joint Wigner function of the final state. As the Wigner function of a state is related to the displaced parity operator, it can be reconstructed by displacing the state and measuring the expectation value of the joint parity gate [12,13,49]. For details on the implementation of the displacement and parity gates, see the review [47].

V. CONCLUSIONS
This article covered the creation of entanglement between bosonic modes in trapped-ion platforms using arithmetic subtractions, a concept first introduced in [28].
First, we presented analytically how arithmetic subtractions can transform a separable state in the local basis into a maximally entangled state over the collective modes. Then, we showed that this scheme can be adapted to transform a state prepared in the collective basis into an entangled state in the same basis. Finally, we presented a numerical simulation of an experimental implementation of this scheme assuming two kinds of noise source: the first due to the coupling to a phasedamping reservoir, and the second to an amplitudedamping reservoir.
In fact, by finding the appropriate beam-splitting operations, this subtraction-based entangling scheme can transform a separable state prepared in any basis (local or collective) into an entangled state in either bases. While we have only considered the case of 4 ions, this scheme can be extended for longer chains. However, since the relations between the local and collective bases for longer chains can only be computed via numerical methods [40], the adaptation of this scheme will need to be done on a case by case basis. On the other hand, the extension of this protocol to the creation of multipartitemultilevel GHZ -like entangled states remains an open problem.
The type of entangled states that can be generated with this subtraction scheme can be used to entangle the internal degree of freedom of different ions with different electronic structures via the red sideband transition. Also, these ions can be trapped in separated wells, and therefore moved to different locations. In addition, the state of one ion can be measured by performing the readout on another ion of a different species. This can beneficial for quantum logic spectroscopy and quantum error corrections [52,53].
In the main text, it was mentioned that the nonzero particle entanglement of the initial state |sym 2n is converted into mode entanglement through the sequence of subtractions. This distinction between mode and particle entanglement arises in the discussion of identical bosonic particles [29], as states appear differently in the first-quantization (particle basis) and secondquantization (occupation number basis) picture. Take, for example, a single-particle basis {|φ 1 , |φ 2 , |φ 3 , |φ 4 }. The state |sym 4 , which describes a system of four particles with one particle in each single-particle state, is written as where we sum over all possible permutations in the particle basis. Then, the mode entanglement of this system is its entanglement in the occupation number basis, and analogously for its particle entanglement. Here, |sym 4 is clearly separable in the occupation number basis, and hence has zero mode entanglement. However, it is highly entangled in the particle basis.
Although the latter form of entanglement is not directly accessible, the notion of particle entanglement is important as it has been shown that nonzero particle entanglement is necessary for extracting mode entanglement from a state using only subtraction operations [30]. The increase in mode entanglement, accompanied by a reduction in the amount of particle entanglement, leads to the statement that one type of entanglement is converted into the other.
However, it is an open problem whether this statement can be read as a quantitative interconversion of resources. To illustrate it, we present an example that takes |sym 4 as initial state and involves only pure states. We quantify the particle entanglement S PE (respectively, mode entanglement S ME ) by calculating the von Neumann entropy maximised over all bipartitions of the state in the particle basis (respectively, occupation number basis). Then obviously for the initial state it holds S PE (|sym 4 ) = log(6), S ME (|sym 4 ) = 0. (A1) We will consider the subtraction sequenceb −θb+θ |sym 4 , whereb +θ andb −θ are defined aŝ The corresponding values of mode and particle entanglement are plotted in Fig. 6. There is indeed a decrease The increase in mode entanglement is accompanied by a decrease in particle entanglement, but the sum SPE + SME is not conserved.
in particle entanglement accompanying the increase in mode entanglement, but the total amount is clearly not conserved.
The scheme being probabilistic, one might think that the probability of success of the transformation should enter quantitative balances. However, this is unlikely to fix the reported discrepancy: for θ → 0, the probability of success tends to 1; nonetheless, the subtraction has significantly reduced S PE while S ME ≈ 0. Perhaps another of the many measures of ME [29] could capture a quantity that is conserved; or perhaps, no quantitative connection should be sought in a transformation between PE and ME. This remains an open question which might be of interest for further theoretical study.
Compared to the target state, this state has the fidelity | ψ f,2n |GHZ 2n | 2 = √ 2 − 1 n + √ 2 + 1 n 2 n (3 n + 1) , and the success probability (3 n + 1)/2 3n−1 . Like before, we can use the red sideband trick to obtain the maximally entangled state. Comparing the state of the system right after the subtraction process, we note that they differ in an extra factor of √ 2 that appears only in the terms where the even modes are oc- Here, g|U rsb,2i π 4 is to be understood as a red sideband transition on the (2i)-th motional mode with θ = π 4 , followed by a post-selection on the ground state of the internal degree of freedom of the ion. This introduces a factor of 1/ √ 2 to the even modes, thus correcting the extra factor of √ 2. The success probability with the red sideband correction is 2 −(2n−1) .
Considering that 2n red sideband transitions are required to prepare the |sym 2n from the motional ground state, this means that a total of 3n first-order gates (red sideband) and 4n second-order gates (beam-splitting and subtraction operations) are required to prepare a |GHZ 2n state for general n.
When considering the total operation time of the general scheme, note that if we coupled lasers to n of the 2n ions in the ion chain, each beam splitter sequence ( n l=1 B † 2l−1,2l , n k=1 B † 2k,2k⊕1 , and n h=1 B 2h,2h⊕1 ) can be performed in parallel. In that case, the time taken for all the beam splitter operations in equation (B8) is 3 × (time taken for one beam splitter operation).