Optimisation of two-dimensional ion trap arrays for quantum simulation

The optimisation of two-dimensional (2D) lattice ion trap geometries for trapped ion quantum simulation is investigated. The geometry is optimised for the highest ratio of ion-ion interaction rate to decoherence rate. To calculate the electric field of such array geometries a numerical simulation based on a"Biot-Savart like law"method is used. In this article we will focus on square, hexagonal and centre rectangular lattices for optimisation. A method for maximising the homogeneity of trapping site properties over an array is presented for arrays of a range of sizes. We show how both the polygon radii and separations scale to optimise the ratio between the interaction and decoherence rate. The optimal polygon radius and separation for a 2D lattice is found to be a function of the ratio between rf voltage and drive frequency applied to the array. We then provide a case study for 171Yb+ ions to show how a two-dimensional quantum simulator array could be designed.


Introduction
Trapped ions possess long lived addressable internal states and can be highly decoupled from their environment. This makes them an important tool in the development of quantum information processing [1,2] and quantum simulation [3][4][5][6][7]. When used for quantum simulation they enable complex spin systems, among others, to be investigated beyond the practical limitations of classical computation. For example, trapped ions have been used for quantum simulations of the evolution of paramagnetic into (anti-)ferromagnetic order in a spin system, [8,9] and frustrated anti-ferromagnetic Ising interactions [10,11]. These first simulations were carried out using one-dimensional trapping arrays and state dependent forces applied using laser beams.
More complex simulations will require ions trapped in two-dimensional (2D) arrays and interaction schemes compatible with these. Advances towards 2D trapping arrays suitable for quantum simulations have been made by trapping ions in a millimetre-scale mechanically fabricated metal mesh [12] and by the successful implementation of microfabrication techniques for ion traps [13]. In addition, interaction schemes based on large oscillating or static magnetic field gradients have been proposed [14,15] and demonstrated in an on-chip microwave gate [16,17]. With the recent advances in the field, ion trap quantum simulations using large scale 2D ion trap lattices are within the reach of current technology. 3 In order to create a 2D array of trapped ions a repeating 2D surface geometry is required. Decoherence due to anomalous heating is a major issue for large scale quantum simulations. As this heating scales approximately as r −4 [18], where r is the ion height above the trapping surface, it is advantageous for ions to be trapped high above the surface. However, when individual surface microtraps are placed together so that their separation is less than around twice the ion height the individual electric fields start to overlap and distort the resulting trapping fields [19]. In extreme cases this can lead to the traps combining to produce a singular trapping zone. To compensate for this the electrode structure has to be altered when operating within this regime [19]. Schmied et al [19] have investigated surface-electrode geometries and developed an algorithm that optimizes geometries to maximize the electric field curvatures of individual trapping sites for arbitrary ion heights and separations. Individual trapping sites shown in [19,20] were optimized using this algorithm leading to non-intuitive electrode patterns which can contain many isolated radio-frequency (rf) and static voltage electrodes. Another proposal [21] working outside this regime uses rf electrodes with controllable rf voltages to lower trap frequencies and decrease ion-ion distances and, therefore, increase interaction strengths. However, this requires the use of multiple independent rf electrodes and individually controllable rf voltages posing an additional experimental complication.
In this paper we present an optimization process for ion trap topologies based on a single rf electrode island, reducing the requirement for buried rf wires and multiple rf electrodes. We focus on the development of an optimum lattice geometry where the ratio of coupling rate to the decoherence rate due to ion heating is maximized and made homogeneous across the lattice. This is achieved by minimizing the secular frequencies of the trapping sites whilst, simultaneously, keeping the trapping depths above a minimum trap depth (for illustrative purposes we use 0.1 eV) to allow for successful operation of the proposed 2D lattice designs. We will concentrate on the optimum lattice topology for hexagonal, square and centred rectangle lattices. An investigation is also carried out on how optimal geometries depend on the overall lattice size, and we discuss the choice of and scaling for experimental parameters such as rf voltage, drive frequency, ion mass and electric field noise density.
Additionally it is possible that 2D ion arrays of this type could also be used for quantum computation. For example, cluster state quantum computation could be carried out in such a system [22,23]. However, additional constraints may also have to be taken into account.

Ion-ion interactions
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ω 2 i , and the Coulomb force, F C = −e 2 /(4π 0 A 2 ), 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]  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] β = e 2 2π 0 mω 2 A 3 , 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.
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 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 The K sim parameter can then be expressed as The secular frequency, ω, of a trapped ion [27] can be expressed as a function of α defined as where V is the amplitude of the rf voltage applied to the trap and is 2π times the drive frequency in Hz, yielding, 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 −4 ω −1 , 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 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).

Two-dimensional ion trap lattice geometry
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. 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].

Simulation of lattices
To determine the electric field produced by a 2D array a method based on the Biot-Savart like law described by Oliveira and Miranda [30] was used. This method calculates the electric field produced by an arbitrarily shaped 2D electrode which is held at a potential V whilst the rest of the plane is held at a potential of zero. The electric field observed at a given point, X , in space due to such an area held at a potential and bounded by a path C is given by [30] where the curve, C, bounds the electrode and x and x are vectors that locate the source point and field point respectively. By calculating the electric field in this manner an assumption is made that there is no gap between the areas held at the potential, V , and the areas held at zero. In microfabricated surface traps, gaps between the electrodes are required and typically range from 3 to 10 µm [13]. If, however, these gaps are small in comparison to the electrode structures they will not alter the trapping fields significantly [20,31,32]. The electric fields of individual electrodes can then be combined to determine potential nils and, therefore, trapping positions in the 2D trap arrays, using the numerical Gauss-Newton algorithm. The secular frequencies, trap depths and ion heights at these positions can then be determined.
To calculate the error of our numerical integration we compared simulations of five wire symmetric surface trap geometries with different central static electrode and rf rail widths in the gapless plane approximation using the method of Oliveira and Miranda [30] to results obtained with analytical equations described by House [31]. In all the geometries simulated the two rf rails were of equal width. Similarly to House, the outer static voltage electrodes were approximated as an infinitely long ground plane although the length of the inner rails were set to 3000 µm instead of infinite. A selection of these simulation results are shown in table 1.
In these results a general error for the ion heights and secular frequencies of less than 2 and 3% respectively was found, which leads to a maximum error in K sim of 10%. For the following simulations it is therefore assumed that the maximum K sim error is 10%. Additionally, numerical simulations of the geometries were carried out using methods described in [33], which indicate similar errors and trends for the ion height and secular frequency as the results obtained with the Biot-Savart like method.

Lattice geometry optimization
In this section we show how the parameters of the lattice geometry (discussed in section 2.2) can be optimized to achieve the highest possible K sim value across the array for a given set of experimental parameters. To do this, we first show how to maximize the homogeneity of individual site properties over an array by varying the distance between the outer polygon in the array to the edge of the rf electrode, g, and show how this scales with lattice size. We will then use the homogeneous arrays to calculate the optimum number of sides, n, a polygon within the array should possess in order to maximize K sim . We then outline a method for the optimization of the polygon radii, R, and separation, A, of an array and show how these vary with increased lattice size and ion mass.

Increasing the homogeneity of K sim across the array
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.
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 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 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 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/ √ N %.   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.

Optimizing the number of polygon sides
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.

Optimization method for polygon separation and radius
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 α 2 and is given by 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.
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 and 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. 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 K sim in equation (9) can also be re-expressed to describe that of an optimized geometry: Equation (19) shows that, for optimal geometries, K sim is proportional to α 3 (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 = (k A /k R ) 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 4 ) [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/α 3 ) 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].

Optimization results and analysis
In this section, optimum polygon separations, A, and radii, R, are obtained using the method outlined above for square, hexagonal and centre rectangular unit cell type lattices. These are shown as functions of lattice sizes and ion masses with the experimental constraint, α, scaled out. Throughout this optimization example, 171 Yb + ions will be used unless otherwise stated. Figure 10 shows how the optimum scaled radius, R/α, and separation, A/α, of polygons vary as a function of lattice size for 171 Yb + ions. As explained in the previous section we have assumed a minimum trap depth of 0.1 eV for illustrative purposes. It can be seen from this figure that as the size of the lattice increases, the optimum polygon radius and separation asymptotically tend towards values representative of an infinitely large lattice. This is expected Figure 11. Graph showing how the optimum K sim /(F 2 α 3 ) varies as a function of the number of sites for optimum lattices with 171 Yb + ions. This is shown for square (square markers), hexagonal (circular markers) and centre rectangular (diamond markers) unit cell lattices. Here F is a state dependent force applied to the ions in the lattice. as once a lattice becomes large enough, the addition of extra lattice sites will represent only a small change in the overall electrode geometry and, therefore, produce a small change in the electric field produced by the geometry. When the lattice is small however, additional lattice sites will represent a larger change in the geometry and will, therefore, cause a larger change in the electric field. Figure 11 shows how the scaled K sim /(F 2 α 3 ) scales as a function of the number of lattice sites, N , using scaled optimum polygon radii, R/α, and separations, A/α. The state dependent force F will be considered in more detail in section 6.2. Due to the dependence of K sim on the geometry, the relationship between K sim /(F 2 α 3 ) and the number of sites is expected to be of similar form to that for optimal polygon radii, R/α, and separation, A/α, with the maximum K sim /(F 2 α 3 ) asymptotically tending towards a value representative of an infinitely large lattice.
Using the data shown in figures 10(a) and (b) and 11 the optimum radii and separation of the polygons as a function of number of sites, N , were found to follow a c + d N −E and f + g N −G relationship, respectively, while the maximum K sim /(F 2 α 3 ) follows a k + l N −Q trend. The values of c, d, E, f , g, G k, l and Q are shown in tables 3-5.
Using the data shown in figures 12(a) and (b) the optimum radii and separation of the polygons as a function of mass were found to follow a o + pm −0.5 and q + sm −0.5 relationship, respectively, with values of o, p, q and s shown in tables 6 and 7. Figure 13 shows how the optimum K sim /(F 2 α 3 ) varies as a function of the mass of the trapped ion, m, for 220 (square Table 4.  type unit cells) and 225 (hexagonal and centre rectangular type unit cells) trapping sites. It is found that the optimum K sim /(F 2 α 3 ) scales as u + vm −0.5 , with the values of u and v shown in table 8.
We note, as shown in figure 12, that as the mass of the ion is increased, the polygon radii and separation will have to be decreased in order to provide trapping regions with a depth of above 0.1 eV for a given α for 220 (square type unit cells) and 225 (hexagonal and centre rectangular type unit cells) trapping sites. It is also clear to see that ions with lighter masses will provide higher K sim /(F 2 α 3 ) values but will require larger lattice geometries compared to heavier ions.  Figure 13. Graph showing how the optimum K sim /(F 2 α 3 ) varies as a function of the ion mass for 220 (square type unit cells (circular markers)) and 225 (hexagonal and centre rectangular type unit cells (square markers and diamond markers respectively)) trapping sites.

Constraints on α
In this section we will discuss the considerations which could limit the value of α. To do this we will show how the power dissipation in a chip trap, the quantum simulation error and the interaction time vary as a function of α. This is important as from this a value of α can be determined for a given experiment, which will be shown in section 7.

Power dissipation in optimized arrays
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] 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] 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 −1 s −1 for the example case of square type lattice with 81 sites using 171 Yb + ions. Considering one particular ion height, r 0 , and substituting for r 0 = k r α 0 and rearranging equation (21) for yields This equation can be re-expressed for V 0 by noting that V 0 = 0 α 0 : Equations (22) and (23) show that, for a given ion mass, m, ion height, r 0 , 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 0 , and unique parameter α 0 . 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 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.

Quantum simulation error
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 Here 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] 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 Equations (19) and (27) show that both the K sim and the error of the simulation, E 0 , 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] 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

20
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 The form of E 0 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 It follows that both the K sim and simulation error E 0 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 171 Yb + will be discussed in section 6.3. Magnetic fields can also be used to provide the state dependent force, F, given by [4] where ω = γ g bi is the position dependent spin resonance frequency with γ g = e/m e the gyromagnetic ratio and i is the x-, yor z-direction. The magnetic field gradient b is assumed to arise from a magnetic field of the form B = B 0 + bî, where B 0 is a constant magnetic field offset. From this, the state dependent force, Fˆi , produced from a magnetic field gradient, bˆi , is found to be 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 Here µ 0 is the permeability of free space, I is the current flowing through the wire and r 2 is the distance squared of the ion from the current carrying wire and is equal to r 2 + a 2 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 The component of this magnetic field gradient in the x-z plane can now be shown to be 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 µ 0 π 2 m 2 e e k 6 r m 2 I 2 cos 2 θ M n + 1 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 α −1 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 α 2 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 ∝ α −1 as it is a function of the inverse of laser sheet cross section, a, which is ∝ α 1 . The magnetic gradient force, on the other hand, is ∝ α −2 as it is a function of the magnetic field gradient b r , which is ∝ α −2 due to r having a linear relationship with α in the geometry considered.

Spontaneous emission
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] 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 171 Yb + ) and g is the single photon Rabi frequency given by Here, I 0 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 It is now possible to describe an additional parameter, L sim , which describes the ratio of interaction rate to the spontaneous emission rate as where It has been shown that the detuning which minimizes the effect of spontaneous emission is approximately 33 THz for 171 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.

Other considerations
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 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), 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 2 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.

Example case study
In this section, an example case will be presented to show how a 2D lattice for the use in quantum simulations can be designed using the work in this paper, whose successful operation is within reach of current technology when using lasers. We then go on to show that while  (23) we can see that there is a unique voltage for a given mass, m, lattice type and stability parameter, q. For this example case we choose a lattice comprised of square type unit cells with nine trapping sites for 171 Yb + ions with a stability parameter q = 0.5. Using these parameters the voltage can now be determined by calculating k r . k r can be found by plotting the ion height of an optimized lattice against α, as shown in figure 8(a), and finding the gradient of this linear relationship. For this example case k r = 98(5) mV −1 s −1 . Using this result and equation (23) we find the unique voltage to be 33(2) V where η geo has been calculated to be 0.14(1).
The α dependent polygon radius and separation can be calculated by producing the corresponding graphs in figure 10, using the method described in section 4.3, for the lattice type and ion to be used. For this example case the optimum radius and separation of the polygons in terms of α are 45(1)α and 174(1)α µm respectively. The next step is to choose a laser detuning from resonance and a maximum acceptable error, E 0 , to be used. For 171 Yb + , as explained in section 6.3, the optimum detuning is approximately 33 THz giving a wavelength of 355 nm, which is therefore used in this example case together with a maximum acceptable error of 0.25. Note that in this example case the array is assumed to be at cryogenic temperatures which correspond to an electric field noise density three orders of magnitude less than at room temperature. We can see from equation (31) that α needs to be minimized in order to keep the quantum simulation error low. The minimum α is determined by the lowest ion height one can easily achieve, which for this example case we choose to be equal to 30 µm. For this case we find α = 0.30(2) V MHz −1 . Having determined α we find the optimum radius and separation to be equal to 14(1) and 52(3) µm respectively. To calculate the laser power required equation (31) should be set to the maximum acceptable error and solved for the laser power, which is found to be 7(1) W assuming the ions to be cooled to n 1 and M = 1 (one site average observable). These conditions provide a coupling rate, J , of 530 ± 110 Hz with β ≈ (2.8 ± 0.2) × 10 −5 . This laser power can, for example, be achieved using a commercially available diode pumped solid state (coherent Paladin range) or fibre (coherent Talisker range) laser system. Table 9 summarizes all parameters for this example case study. Figure 14 shows the effect a change of α and laser power, P, has on the quantum simulation error (solid curves), K sim (dashed curves) and L sim (dotted lines). The cross corresponds to the 2D lattice designed in this example case which represents the optimum case in terms of achieving the highest K sim and L sim for a maximum quantum simulation error of 0.25. We note that the main limitations in achieving lower quantum simulation errors stem from the lowest achievable ion height (lowest α) and magnitude of electric field noise density. Figure 14 also shows that higher values of K sim and L sim can be achieved with any given geometry (given by the value of α) by simply increasing the power of the laser. However, this can only be achieved at the expense of higher quantum simulation errors.
State dependent forces can also be created using magnetic field gradients as described in section 6.2. Figure 15 shows the K sim (solid curves) and the quantum simulation error (dashed curves) as a function of the magnetic field gradient, b, and α for traps operated at cryogenic temperature. Here we use 171 Yb + ions in a 3 × 3 square unit cell array. As described in section 6.2, E 0 ∝ α −1 indicating that a large α is advantageous. The limit on the maximum α is dependent on the maximum current one can apply to the geometry. Using equation (35) it is possible to calculate the current required to create a desired magnetic field gradient. In order to illustrate the magnitude of currents required we assume k a = k r , which will result in an angle θ = 45 • (refer to section 6.2 for more information). We also set α = 0.30(2) V MHz −1 , determined by the lowest achievable ion height which, for illustration purposes, we have set to 30 µm. The reason for choosing the minimum α value can be seen when considering equation (35) which clearly shows that, for a given magnetic field gradient, I ∝ α 2 . In the magnetic field gradient case, the chosen α sets K sim and E 0 . For this case, again, we choose M = 1, n 1, K sim = 2 and E 0 = 0.01 which requires a magnetic field gradient of 30 500 ± 3500 T m −1 and is indicated by the cross in figure 15. This is achievable with a current of approximately 1200 A yielding a coupling rate, J , of 240 Hz and a β = (2.8 ± 0.2) × 10 −8 . From this simple example case one can conclude that using magnetic field gradients to provide state dependent forces for use in quantum simulations using the methods and trap designs shown in this work is quite challenging. However, geometries trapping ions in chains allow for sizeable magnetic gradient induced couplings [39] and, so, a detailed investigation into optimizing the wires used for producing magnetic field gradients in the geometries discussed in this work could improve results.

Conclusion
2D arrays of surface ion traps have the potential to provide a technology with which quantum simulations can be performed. In order for ion traps to be used successfully for this purpose a greater understanding is required of how the various geometry parameters affect the ions trapped above them. Throughout this work square, hexagonal and centre rectangular unit cell arrays of microtraps have been modelled in the gap-less plane approximation using the Biot-Savart like law in electrostatics [30]. Decoherence due to motional heating [18] was then compared to the ion-ion interaction [24] to provide a ratio used to describe how much faster an ion-ion interaction occurs in comparison to the motional decoherence, K sim . This work investigates how various parameters in the array can be adjusted in order to optimize the device's ability to perform quantum simulations and shows how the interactions can be made as homogeneous as possible over the device. It has been shown how the homogeneity of the K sim across an array can be altered by varying the distance from the outer polygon to the edge of the rf electrode. The distance required to maximize the K sim homogeneity is also shown to vary as a function of the total size of the lattice. The number of polygon sides, n, required to maximize K sim has also been found.
We have shown that the K sim of a given lattice geometry can be maximized by reducing the value of α. However, as α reduces so does the trap depth. This results in the conclusion that the maximum K sim of a geometry can be achieved by reducing the value of α until the trap depth reaches a reasonable minimum value.
Using this information, optimal geometries as a function of α are presented. This has been achieved by finding the relationships of polygon separation and radius with α for optimal geometries. It was found that, for these optimal geometries, K sim scales as α 3 . The individual polygon separation and radius were found to possess a linear relationship with α and, therefore, larger geometries have been found to produce larger values of K sim . Therefore, the optimal lattice geometry is dependent solely on the value of α for a given ion mass and number of trapping sites in the array.
We presented a case study for determining an optimum geometry consisting of nine trapping sites arranged into square type unit cells for 171 Yb + ions. We showed how the value of α can be chosen (and, therefore, the geometry dimensions) by taking into account the laser power (or static magnetic field gradient) required to produce a state dependent force acting on the ions, the K sim and the error on the simulation. From this it has been found that to carry out quantum simulations with reasonable K sim and error values it is preferable to use traps held at cryogenic temperatures as this reduces decoherence due to heating effects on the ions. Other methods known to significantly decrease the anomalous heating rate include pulsed laser electrode cleaning [40] and argon-ion beam electrode cleaning [41].
The scaling of anomalous heating has not yet been fully understood and is thought to possess a dependence on the ion height, r , of between r −4 and r −2 . In this work we have used r −4 , however, if the anomalous heating is found to possess a relationship with ion height nearer r −2 then the equations in this work can be adjusted and this is discussed in more detail in appendix B.
The relationships between lattice size and α with the polygon radii and separation obtained using the method described in this work will allow for the construction of 2D surface trap lattice arrays with high ratios of ion-ion interaction rates to decoherence rates, providing a system which could be used to perform quantum simulations.
The simulation error in equation (27) can be adjusted to take into account the heating of the ions during a simulation. If this occurs the mean radial phonon number n will become time dependent, n(t). It, therefore, follows that the error will also become time dependent and is given by The time dependent mean radial phonon number will be a function of the heating rateṅ, the interaction time, T J , and the initial mean phonon number, n 0 , and is given by n(T J ) = n 0 +ṅT J (A.2) yielding an error given by E 0 (T J ) ≈ F 2 M(n 0 +ṅT J + 1 2 ) 2hmω 3 .
(A.4) Equation (A.4) shows that as K sim increases the error will tend towards that in equation (27) in section 6.2 for the case ofn 1. This is because higher values of K sim mean that less heating takes place during an interaction until the mean radial phonon number can be approximated as constant.

Appendix B
In the work presented we have used a motional heating rate,ṅ, which has an r −4 scaling, where r is the ion height. However the scaling of this anomalous heating has not yet been fully understood and so it is conceivable that a different scaling will be required to describe the effect. For this reason it is the aim of this appendix to outline the steps and main expressions required to obtain an optimized 2D ion trap array with an anomalous heating rate which possesses an arbitrary scaling with the ion height r −x .
For the case of optimizing the homogeneity of K sim across an array the results shown in section 4.1 will hold for any scaling of the heating rate with ion height. This is because this optimization aims to give each trapping site in the array the same properties and this is independent of the heating rate. The same applies to the optimization of the number of polygon sides described in section 4.2.
Equation (19) describing K sim for an optimized geometry can be altered to take into account a different scaling of the heating rate with ion height. For an arbitrary scaling r −x this equation can be expressed as It can be seen in equation (B.1) that K sim is proportional to α (x−1) . The values of k r , k A and k R (described in section 4.3) are independent of the value of x and, so, the optimum geometry for a given α will be the same regardless of the scaling of the heating rate with ion height. Using this procedure optimum geometries can be computed for an arbitrary scaling of the anomalous heating with the ion height, r −x .