Optimization of two-dimensional ion trap arrays for quantum simulation

2D lattices of ions can be used as a quantum simulator for many body spin-1/2 systems [24]. Forces such as the trapping potential, F_T = −mω²_i, and the Coulomb force, F_C = −e²/(4π₀A²), between the individual ions determine the equilibrium position. Laser beams or magnetic field gradients can be used to impart an additional force to the ions which displace the ion(s) depending on its internal state. This displacement leads to a change in the Coulomb force and thereby displaces the neighbouring ion(s) dependent on their own internal state. This state dependent coupling is given by [24]

$$\begin{equation} J = \frac{\beta F^{2}}{m\omega^{2}}, \end{equation} \tag{ 1 }$$

where m is the mass of an ion, F is the magnitude of the state dependent force applied to each ion and ω is the trap's secular frequency. We will consider how this force is applied later in the article. β is the ratio of the change in the Coulomb force to the change in the restoring force due to the displacement of the ions caused by the state dependent force and is given by [24]

$$\begin{equation} \beta = \frac{e^{2}}{2\pi\epsilon_{0}m\omega^{2}A^{3}}, \end{equation}\noindent \tag{ 2 }$$

where A is the ion–ion spacing. There are two cases to consider, when β > 1 the change in the Coulomb force, δF_C, due to the displacement of an ion is dominant over the change in the restoring force, δF_T. This results in an interaction over a large number of trapping sites. When β < 1, the opposite is true resulting in an interaction which decays rapidly across the array. Trapping ions in a 2D array of microtraps makes it possible to satisfy the condition that β < 1 allowing systems with short range interactions to be simulated. An illustration of this system is shown in figure 1.

Zoom In Zoom Out Reset image size

It is important to consider sources of decoherence when designing a 2D ion trap array. The internal state of an ion can remain coherent for tens of seconds [25, 26]. However, motional decoherence due to anomalous heating of ions will be an important factor during quantum operations within small scale ion traps as the implementation of spin dependent couplings involves the use of motional states of the ion. The coupling, J, will be observable if the coupling time, T_J = 1/J, is less than the motional decoherence time in the system and, therefore, the ratio of these two times is an important parameter of the system and is given by

$$\begin{equation} K_{\mathrm{sim}}=\frac{T_{\dot{n}}}{T_{J}}, \end{equation} \tag{ 3 }$$

where [18]

$$\begin{equation} T_{\dot{n}}=\frac{4m\omega\hbar}{e^{2}S_{\mathrm{E}}(\omega)}. \end{equation} \tag{ 4 }$$

Here S_E(ω) is the electric field noise density [13, 18]. In order for an interaction to occur on faster time-scales than the decoherence in the system, we require K_sim > 1 and it is the aim of the optimization process presented in this work to optimize the geometry in order to maximize this parameter. To acquire an understanding of how a geometry can affect K_sim it is necessary to determine its form with respect to the geometry variables. The form of T_J can be found by substituting equation (2) into (1) and is given by

$$\begin{equation} T_J=\frac{2\pi\epsilon_{0}m^{2}\omega^{4}A^{3}\hbar}{e^{2}F^2}. \end{equation} \tag{ 5 }$$

The K_sim parameter can then be expressed as

$$\begin{equation} K_{\mathrm{sim}}=\frac{2F^2}{S_{\mathrm{E}}(\omega)\pi\epsilon_{0}m\omega^{3}A^{3}}. \end{equation} \tag{ 6 }$$

The secular frequency, ω, of a trapped ion [27] can be expressed as a function of α defined as

$$\begin{equation} \alpha=\frac{V}{\Omega}, \end{equation} \tag{ 7 }$$

where V is the amplitude of the rf voltage applied to the trap and Ω is 2π times the drive frequency in Hz, yielding,

$$\begin{equation} \omega=\frac{eV\eta_{\mathrm{geo}}}{\sqrt{2}m\Omega r^2}=\frac{e\alpha\eta_{\mathrm{geo}}}{\sqrt{2}m r^2}, \end{equation} \tag{ 8 }$$

where r is the height of an ion above the surface, e is the charge of an electron and η_geo is an efficiency factor which can range between 0 and 1 depending on the form of the geometry [27].

The secular frequency given in equation (8) can then be used along with S_E(ω) = Ξr⁻⁴ω⁻¹, where Ξ is a coefficient that can be experimentally obtained and depends on the temperature and surface of the trap electrodes (see [13] for a listing) to re-express equation (6) in the form

$$\begin{equation} K_{\mathrm{sim}}=\frac{4F^2mr^8}{\Xi\alpha^2\eta_{\mathrm{geo}}^2\pi\epsilon_0A^3}. \end{equation} \tag{ 9 }$$

To further understand how the geometry effects the K_sim value we will now introduce the parameters of a lattice geometry and relate them to equation (9).

A lattice is a regular tiling of a space by a primitive unit cell. Previous works [12, 21] concentrate on lattices created from square unit cells. In total there exist five types of cell which can be used to form a 2D lattice: centred rectangular, hexagonal and square as shown in figure 2, and rectangular and oblique [28]. The rectangular and oblique structures are not considered in this work due to their non-uniform ion–ion distances.

Zoom In Zoom Out Reset image size

Figure 2 shows the polygon–polygon separation which is equal to the ion–ion distance, A, in equation (9). The polygon radius, R, along with the separation, will determine the height above the surface at which the ion is trapped, r, with larger polygon radii yielding higher ion heights. Another variable to be considered is the gap between the outer polygon in the array to the edge of the rf electrode, g. This can be used to alter the homogeneity of the individual trapping sites within the array. In general a non-homogeneous system results in spin dependent coupling rates which are a function of the lattice site, posing a significant problem for the scalability of such an array [29].

This is achieved by ensuring homogeneous secular frequencies, ion heights and trap depths across all the array sites. As shown in figure 3 the K_sim of trapping sites in an array can be altered to approach a common value if the distance, g, between the edge of the outer polygon and the edge of the rf electrode containing the polygon array is adjusted. As the value of g is increased, the K_sim value of the sites towards the centre drops, however, the outer sites K_sim value rises, resulting in the properties converging towards each other. If the distance g is increased further beyond the point at which maximum homogeneity occurs the outer site properties drift away from those of the central sites and, therefore, decrease the homogeneity of the array.

Zoom In Zoom Out Reset image size

To provide a value of g which is universal for all lattice sites, its value is given in units of lattice side length, L. The lattice side length is determined by the polygon separation, A, and radius, R, and can be expressed as

$$\begin{equation} L=(M-1)A+2R, \end{equation} \tag{ 11 }$$

where M is the number of lattice sites along one side of an array.

In order to quantify the array's homogeneity, H is defined as the average deviation of K_sim of each lattice site from the K_sim of the central site and is given by

$$\begin{equation} H = \frac{1}{N} \sum_{n=1}^{N} \left|1-\frac{K_{\mathrm{sim}_{n}}}{K_{\mathrm{sim}_{\mathrm{centre}}}}\right|, \end{equation} \tag{ 12 }$$

where N is the total number of trapping sites in the lattice.

Figure 4 shows H for a 5 × 5 square type unit cell lattice for 0 < g/L < 1.5. The maximum homogeneity, and thus the optimum g/L, is found when H is minimized. The error associated with H is given by

$$\begin{equation} \sigma_{H} = \frac{\sqrt{N(\sigma_{K\mathrm{sim}})^{2}+N(\sigma_{K\mathrm{sim}_{\mathrm{centre}}})^{2}}}{N}, \end{equation} \tag{ 13 }$$

where the error on all K_sim values is 10%, as shown in section 3. This yields an overall percentage error on H of 0.13/$\sqrt {N}\%$.

Zoom In Zoom Out Reset image size

Figure 5 shows the optimum g/L for hexagonal, central rectangular and square unit cell lattices of different sizes. The curves are found to be described by an equation of the form g/L = a + bN^−B with a, b and B values for different lattice types shown in table 2. For large lattices, g/L is independent of N, as trapping sites close to the edge of the lattice are influenced only by the electric field created by that edge. In small lattices, however, the optimum g/L increases as the effect of the electric field from the opposite edge of the lattice increases.

Zoom In Zoom Out Reset image size

We now investigate the optimum number of polygon sides, n, providing the highest K_sim value on the central trapping site (located above the central polygon) of a lattice.

To ensure the results are universal for all lattice geometries, g/L is set to the value which maximizes the homogeneity of each array, and all other parameters are scaled by normalizing K_sim to that of an identical geometry with polygons of 100 sides. This also allows comparison between the different types of lattices.

Figure 6 shows the scaled K_sim for the central site as the number of polygon sides is varied. As the number of sides is increased the value of K_sim approaches an asymptote, shown by the upper dashed line. This indicates that the best geometry will be made from circular electrodes. However, simulation times grow with increasing polygon side number and so it is advantageous to reduce this number to a minimum. It is shown that 95% of the maximum achievable K_sim (indicated by the lower dashed line) can be achieved with around ⩾25–30 sides in the polygons.

Zoom In Zoom Out Reset image size

In this section we now maximize the K_sim of any arbitrary geometry. We will then go on and determine optimum geometries and show how they scale and, ultimately, are determined by α = V/Ω. We use the value of g determined in section 4.1 in order to provide the maximum K_sim homogeneity across the lattice, and set the number of sides to 25 as this provides a good approximation to the optimum circular geometry while keeping the simulation time at a minimum, as shown in section 4.2.

When considering any fixed arbitrary geometry, equation (9) shows that K_sim can be maximized by reducing the value of α. However, the minimum achievable α is limited by the lowest usable trap depth, as the trap depth is proportional to α² and is given by

$$\begin{equation} T_{\mathrm{D}}=\frac{\zeta e^2\alpha^2}{\pi^2m}, \end{equation}\noindent \tag{ 14 }$$

where e is the charge of an ion and ζ is a geometrical factor which is a function of A and R [31].

We will now focus on finding optimal geometries which we define as geometries that yield the highest value of K_sim for a given value of α. This will be carried out by fixing the trap depth at a reasonable minimum value (we use 0.1 eV for illustration purposes) which, as discussed above, provides the maximum K_sim for a given geometry, and then investigating the dependence of polygon radius, separation and ion height with α. To determine these dependences a K_sim contour plot is made by calculating the K_sim over a range of polygon separation, A, and radius, R, with a resolution of 1 μm. The range of polygon radius and separation used should not create traps with inter-well barriers of less than the minimum trap depth value, and to ensure this the polygon radius was kept to less than a third of the polygon separation. For each combination of polygon separation and radius a value of α is found which yields the minimum trap depth and, thus, maximizes the K_sim of the particular geometry. By following this method one can obtain the α required to achieve the minimum trap depth, the ion height, r, and K_sim for each geometry. From the resulting data the polygon separation and radius which yields the highest K_sim for a given α (the optimum geometry) can then be found. A graphical example of such data is shown in figure 7.

Zoom In Zoom Out Reset image size

It can be seen from this method and the example data in figure 7 that the highest K_sim will be achieved with an infinite value of α, R and A. However, other effects may limit the magnitude of α. In order to determine a limit on α it is, therefore, necessary to describe the various array and trapping field dependent properties (such as secular frequency, ion height and K_sim) in terms of α.

By plotting these optimum parameters (polygon radius, separation and ion height) as a function of α, as shown in figures 8(a)–(c), linear relationships of the form

$$\begin{equation} r=k_r\alpha, \end{equation} \tag{ 15 }$$

$$\begin{equation} A=k_{A}\alpha \end{equation} \tag{ 16 }$$

and

$$\begin{equation} R=k_{R}\alpha \end{equation}\noindent \tag{ 17 }$$

are found for the optimal geometries. The values of k_r, k_A and k_R are dependent on the number of trapping sites in an array, as shown in figures 9(a)–(c) respectively, for lattices made from square type unit cells of polygons.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is important to stress that equations (15)–(17) are only valid in the case of optimal geometries, which depend solely on α. With this in mind, it is now possible to re-express the secular frequency in equation (8) to describe the secular frequency of an optimized geometry as

$$\begin{equation} \omega=\frac{e\eta_{\mathrm{geo}}}{\sqrt{2}k_{r}^{2}m\alpha}. \end{equation} \tag{ 18 }$$

K_sim in equation (9) can also be re-expressed to describe that of an optimized geometry:

$$\begin{equation} K_{\mathrm{sim}}=\frac{4F^2mk_{r}^{8}\alpha^3}{\Xi\eta_{\mathrm{geo}}^2\pi\epsilon_{0}k_{A}^{3}}. \end{equation} \tag{ 19 }$$

Equation (19) shows that, for optimal geometries, K_sim is proportional to α³ (as the value of α determines the electrode dimensions to be used) and so, to produce an array with a high K_sim for a given number of lattice sites (as k_r and k_A are functions of the number of trapping sites), a large value of α is preferable. Equations (16) and (17) show that the optimum geometry size is proportional to α. It, therefore, follows that larger lattice geometries will produce larger values of K_sim. This effect is illustrated in figure 7 which shows the K_sim as a function of polygon radius and separation. For optimized lattices, the optimal radius and separation will fall on a line described by $A=\left (k_A/k_R\right )R$, with higher values of α required for higher values of separation and radius as shown in equations (16) and (17) respectively.

The heating rate in ion traps has a strong dependence on the ion height (∝1/r⁴) [18]. A large K_sim is achieved with large values of α, resulting in large ion heights, as shown in equation (15). It, therefore, can be concluded that a different scaling of the heating rate (for example in cryogenic systems) does not change the optimization process and optimal geometries.

It has now been shown that the optimal geometry for a given minimum trap depth and ion mass is determined solely by the value of α.

Optimal geometries and their K_sim values (in units of 1/α³) can now be found by creating contour plots (such as shown in figure 7) for different values of lattice size, N, lattice unit cell type and ion mass, m. The error on the K_sim value, calculated from equation (9), was determined to be ±10% by comparing the secular frequency and ion heights obtained using the program with those predicted by House's analytical solutions for a five wire surface trap geometry [31].

The power dissipation of an ion trap chip is determined by the voltage, V , frequency, Ω, as well as the capacitance and resistance of the trap itself. This may, for a given capacitance and resistance, affect the value of α (the ratio between the rf voltage and drive frequency) which can be used. As the value of α is used to determine the optimum polygon radii and separation of a geometry, as shown in figure 10, for a given number of sites, N, and stability parameter, q, it is important to know how the power dissipation varies as a function of α.

The power dissipation of a chip is approximately given by [13]

$$\begin{equation} P_{\mathrm{D}}\approx\textstyle\frac{1}{2}V^{2}C^{2}R\Omega^{2}, \end{equation} \tag{ 20 }$$

where C and R are the capacitance and resistance of the chip. It is not possible to apply any combination of V and Ω to a geometry as they must be chosen so that the ion is stably trapped with a stability parameter given by [12, 27]

$$\begin{equation} q=\frac{2e\eta_{\mathrm{geo}} V}{mr^2\Omega^2} =\frac{2e\eta_{\mathrm{geo}}\alpha}{mr^2\Omega} \end{equation} \tag{ 21 }$$

between 0 and 0.9, where e is the charge of an electron.

The ion height, r, of ions trapped in the optimized lattices shown in this work have been found to be linearly proportional to α. This relationship is shown in figure 8(a) with the constant of proportionality, k_r, found to be 63(4) mV⁻¹ s⁻¹ for the example case of square type lattice with 81 sites using ¹⁷¹Yb⁺ ions. Considering one particular ion height, r₀, and substituting for r₀ = k_rα₀ and rearranging equation (21) for Ω yields

$$\begin{equation} \Omega_0=\frac{2e\eta_{\mathrm{geo}}}{mk^2_r\alpha_0 q}. \end{equation} \tag{ 22 }$$

This equation can be re-expressed for V₀ by noting that V₀ = Ω₀α₀:

$$\begin{equation} V_0=\frac{2e\eta_{\mathrm{geo}}}{mk^{2}_{r}q}. \end{equation} \tag{ 23 }$$

Equations (22) and (23) show that, for a given ion mass, m, ion height, r₀, stability parameter, q, and number of trapping sites in the array (as k_r is a function of the number of trapping sites), there is one unique voltage, V₀, and unique parameter α₀. This means that a given ion height (and, therefore, a chosen value of α) determines both the voltage and drive frequency to be applied to the trap.

To express the power dissipation, P_D, in terms of α equations (22) and (23) can be substituted into equation (21) giving

$$\begin{equation} P_{\mathrm{D}}=\frac{8\eta_{\rm geo}^4e^4C^2R}{k_r^8m^4q^4\alpha^2}. \end{equation} \tag{ 24 }$$

Equation (24) shows that as α is increased, the power dissipated is reduced. This means that the power dissipation is low for high values of α and, as high values of α provide high values of K_sim (see figures 11 and 13), power dissipation will not impact on producing high values of K_sim in optimized geometries. However, a low value of α will result in a high power dissipation in the chip and, so, the maximum allowable power dissipation in a chip will determine the lowest α which can be applied to a geometry.

An upper limit on the value of α can be obtained from an estimation of the error of a quantum simulation using the method described in [34], where the error for the Ising model is given by

$$\begin{equation} E_0\approx\frac{1}{2}\eta^2\sum_j(2\overline{n}+1)\left\langle\left[\left[O(t),\sigma_j^z(t)\right],\sigma_j^z(t)\right]\right\rangle. \end{equation} \tag{ 25 }$$

Here $\overline {n}$ is the mean radial phonon number of the ions, O is the observable of the quantum simulation and η is a parameter which characterizes phonon displacement caused by the state dependent force and is given by [34]

$$\begin{equation} \eta=\frac{F\sqrt{\hbar/(2m\omega)}}{\hbar\omega}, \end{equation} \tag{ 26 }$$

where m is the mass of a trapped ion and ω/2π is its secular frequency.

If O is an M-site observable then there exist M non-vanishing commutators (for example a two-site correlation function (M = 2) or a spin average (M = 1)) and so the error on the simulation will not be dependent on the number of ions, N, in the array [34]. The error in equation (25) can now be re-written as

$$\begin{equation} E_0\approx\frac{F^{2}M(\overline{n}+\frac{1}{2})}{2\hbar m\omega^{3}}. \end{equation} \tag{ 27 }$$

Equations (19) and (27) show that both the K_sim and the error of the simulation, E₀, are proportional to the square of the state dependent force, F. It follows that the way in which this force is applied to the ions will determine the dependence of the simulation error on α. A discussion on the possible effects of the heating rate on the error can be found in appendix A. In this work we will consider applying this state dependent force via laser beams and magnetic field gradients.

To calculate the laser power required to achieve a force, F, it is assumed, for illustrative purposes, that the laser beam is focused to a sheet given by 25 μm multiplied by the width of the array. The laser intensity required to provide a state dependent force, F, can be provided by a laser beam of power, P. If the output power of the laser used is assumed to be constant, then the force applied to the ions will be dependent on α. This is because the lattice size will increase with increasing α and, therefore, so will the spatial area of the beam required to impart a force onto the ions. It is, therefore, required to express this force as a function of α. The intensity of a beam required to provide a state dependent force, F, is given by [35]

$$\begin{equation} I_0=\frac{3F\Delta\lambda I_{\mathrm{sat}}}{2\pi\hbar\gamma^2}, \end{equation} \tag{ 28 }$$

where Δ is the detuning of the laser from resonance, λ is the wavelength of the laser, I_sat is the saturation intensity of the ion and γ is 2π times the transition linewidth. The power of a laser beam is given by

$$\begin{equation} P=I_0a, \end{equation} \tag{ 29 }$$

where a is the spatial area to which the beam is focused. In this work the beam is assumed to be focused to form a light sheet across the array with an area given by a = (n_s − 1)AW = (n_s − 1)k_AαW, where n_s is the number of trapping sites (or polygons) along one side of the array and W is the width of the light sheet. By using equations (28) and (29) the force applied to the ions by a laser power, P, can be expressed as

$$\begin{equation} F=\frac{2\pi\hbar P\gamma^2}{3a\Delta\lambda I_{\mathrm{sat}}}. \end{equation}\noindent \tag{ 30 }$$

The form of E₀ for the case of lasers applying the state dependent force can now be found by using equations (18), (30) and the general error equation (27) yielding

$$\begin{equation} E_{0 \mathrm{laser}}=\frac{4\sqrt{2}}{9}\frac{\pi^2\hbar}{e^3}\frac{M k_r^6\gamma^4m^2P^2\alpha(\overline{n}+\frac{1}{2})}{\eta_{\mathrm{geo}}^3(n_{\mathrm{s}}-1)^2k_A^2W^2\Delta^2\lambda^2I_{\mathrm{sat}}^2}. \end{equation}\noindent \tag{ 31 }$$

It follows that both the K_sim and simulation error E₀ will decrease with increasing laser detuning, Δ. Therefore, the optimum detuning, for a given laser power and α, corresponds to the lowest detuning which provides the required K_sim. The optimum detuning to achieve the lowest simulation error for ¹⁷¹Yb⁺ will be discussed in section 6.3.

Magnetic fields can also be used to provide the state dependent force, F, given by [4]

$$\begin{equation} F_{\hat{i}}=\left(\frac{\hbar}{2}\right)\partial_{i}\omega\langle \sigma^{(\hat{i})}\rangle, \end{equation}\noindent \tag{ 32 }$$

where ω = γ_gbi is the position dependent spin resonance frequency with γ_g = e/m_e the gyromagnetic ratio and i is the x-, y- or z-direction. The magnetic field gradient b is assumed to arise from a magnetic field of the form $\overline {B}=\overline {B_0}+b\hat {i}$, where B₀ is a constant magnetic field offset. From this, the state dependent force, $F_{\hat {i}}$, produced from a magnetic field gradient, $b_{\hat {i}}$, is found to be

$$\begin{equation} F_{\hat{i}}\approx\frac{\hbar eb_{\hat{i}}}{2m_{e}}, \end{equation}\noindent \tag{ 33 }$$

where m_e and e are the mass and charge of an electron respectively. If one assumes the magnetic field is created by a current carrying wire located on the surface of a polygon, at a distance a from the centre of the polygon and making an angle of 45° with respect to the x-axis, then the magnetic field gradient will be of the form

$$\begin{equation} b_{r^{\prime}}=\frac{\mu_0I}{2\pi r^{\prime 2}}. \end{equation}\noindent \tag{ 34 }$$

Here μ₀ is the permeability of free space, I is the current flowing through the wire and r'² is the distance squared of the ion from the current carrying wire and is equal to r² + a² where r is the ion height. We assume, for simplicity, that the distance a scales linearly with α with a constant of proportionality of k_a in order to keep the angle of r' to the x–z plane, θ, independent of α. As the ion height is known to scale linearly with α, from equation (15), it is possible to express the magnetic field gradient along r' as

$$\begin{equation} b_{r^{\prime}}=\frac{\mu_0I}{2\pi\alpha^2\left(k_r^2+k_a^2\right)}. \end{equation} \tag{ 35 }$$

The component of this magnetic field gradient in the x–z plane can now be shown to be

$$\begin{equation} b_{x,z}=\frac{\mu_0I\cos\theta}{4\pi\alpha^2\left(k_r^2+k_a^2\right)}. \end{equation} \tag{ 36 }$$

The form of K_sim for the case of magnetic field gradients applying the state dependent force can be found by using equations (18), (33), (36) and the general error equation (27) yielding

$$\begin{equation} E_{0 \mathrm{mag}}=\frac{\sqrt{2}}{64}\frac{\hbar\mu_0}{\pi^2m_e^2e}\frac{k_r^6m^2I^2\cos^2\theta M\left(\overline{n}+\frac{1}{2}\right)}{\eta_{\mathrm{geo}}^3\left(k_r^2+k_a^2\right)^2\alpha}. \end{equation}\noindent \tag{ 37 }$$

Equations (31) and (37) show that the quantum simulation error is proportional to α for a state dependent force created by a laser beam and proportional to α⁻¹ for a magnetic field gradient created by current carrying wires. For the case of laser beams the α scaling implies that as α is increased (yielding larger K_sim values and geometries as shown in section 4.3) the quantum simulation error will rise and, therefore, provide an upper limit on the value of α. For the magnetic field gradient case the upper limit on α comes from the current creating the gradient. While the α scaling for the quantum simulation error using magnetic field gradients implies that a larger α is advantageous, the current required to achieve a given magnetic field gradient scales as α² as can be deduced from equation (35). The maximum current that one can apply to the lattice therefore provides an upper limit for α.

It is also interesting at this point to note the different scaling with α of the laser and magnetic field gradient forces given in equations (30) and (32). The laser force can be seen to be ∝α⁻¹ as it is a function of the inverse of laser sheet cross section, a, which is ∝α¹. The magnetic gradient force, on the other hand, is ∝α⁻² as it is a function of the magnetic field gradient b_r', which is ∝α⁻² due to r' having a linear relationship with α in the geometry considered.

When applying the state dependent force to the ions using a laser beam additional decoherence will occur via spontaneous emission. The spontaneous emission rate is given by [36, 37]

$$\begin{equation} S=\frac{\gamma g^2}{6}\left(\frac{1}{\Delta^2}+\frac{2}{\left(\Delta_{\mathrm{fs}}-\Delta\right)^2}\right), \end{equation}\noindent \tag{ 38 }$$

where γ is 2π times the linewidth in Hz, Δ is the laser detuning from resonance, Δ_fs is the fine structure splitting of the ion (approximately 100 THz for ¹⁷¹Yb⁺) and g is the single photon Rabi frequency given by

$$\begin{equation} g=\gamma\sqrt{\frac{I_0}{2I_{\mathrm{sat}}}}. \end{equation} \tag{ 39 }$$

Here, I₀ is the laser intensity and I_sat is the saturation intensity of the ion. It is possible to express the single photon Rabi frequency in terms of the laser power, P, by using equation (29) giving

$$\begin{equation} g=\gamma\sqrt{\frac{P}{2(n_{\mathrm{s}}-1)k_A\alpha WI_{\mathrm{sat}}}}. \end{equation} \tag{ 40 }$$

It is now possible to describe an additional parameter, L_sim, which describes the ratio of interaction rate to the spontaneous emission rate as

$$\begin{equation} L_{\mathrm{sim}}=\frac{T_S}{T_J}, \end{equation} \tag{ 41 }$$

where

$$\begin{equation} T_{S}=\frac{1}{S}. \end{equation}\noindent \tag{ 42 }$$

It has been shown that the detuning which minimizes the effect of spontaneous emission is approximately 33 THz for ¹⁷¹Yb⁺ [38]. It, therefore, follows that the value of L_sim will be maximized with this detuning. This additional parameter is analogous to the parameter K_sim in equation (3) and is also required to be greater than unity, just like the original K_sim, when considering a state dependent force created using laser light. If magnetic field gradients are to be used to apply the state dependent force then L_sim is no longer relevant.

It is important to note here that an increase in α will increase the time taken for ion–ion interactions to take place in optimized lattice structures, as we will show in equation (44). Equation (5) gives an expression for the time taken for an ion–ion interaction to occur in any given fixed lattice structure. This can be expressed for optimized lattice structures by including the expressions for the ion–ion separation (polygon separation, A) and secular frequency, ω, from equations (16) and (18) respectively, to yield

$$\begin{equation} T_J=\frac{e^2\pi\epsilon_0\hbar}{2}\frac{\eta_{\mathrm{geo}}^4k_A^3}{k_r^8F^2m^2\alpha}. \end{equation}\noindent \tag{ 43 }$$

An expression to give the interaction time in optimized lattices as a function of α can now be arrived at by using equations (43) and (30),

$$\begin{equation} T_J=\frac{9}{8}\frac{e^2\epsilon_0}{\pi\hbar}\frac{\eta_{\mathrm{geo}}^4k_A^3((n_s-1)k_AW)^2\Delta^2\lambda^2I_{\mathrm{sat}}^2\alpha}{k_r^8\gamma^4P^2m^2}. \end{equation}\noindent \tag{ 44 }$$

Equation (44) clearly shows that as α is increased the time taken for an ion–ion interaction will increase and, so, it may be preferable to limit the magnitude of α after taking into account the effects on K_sim. A similar equation can also be derived for use with magnetic field gradients. It is also important to note here that increasing the laser power, P, will increase the value of L_sim as the spontaneous emission rate is proportional to P whereas the coupling rate is proportional to P² as shown in equations (38) and (44) respectively.

With the use of the equations derived in this section, optimal geometries can be calculated given certain experimental parameters and will be described in the following section.