{'fraction_non_alphanumeric': 0.03419161990665478, 'fraction_numerical': 0.009142237838026816, 'mean_word_length': 4.201673031460266, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 22, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Perturbative approach to quantum tunneling and quantum interferences in spin clusters

Collective tunneling is a ubiquitous phenomenon in finite-size spin clusters that shows up in systems as diverse as molecular magnets or spin clusters adsorbed at surfaces. The problem we explore is to understand how small flipping terms can cooperate to flip a large spin to the opposite direction or a cluster of interacting elementary Ising spins into the time-reversed state. These high order processes will involve at least two channels, a single spin-flip channel due to a transverse field and a two-spin flip channel due to exchange or due to single-ion anisotropies. In the present paper, we show that high-order perturbation theory can be formulated and evaluated with the help of simple recurrence relations, leading to a compact theory of tunnelling in macroscopic spins, in one-dimensional clusters, as well as in small higher-dimensional clusters. This is demonstrated explicitly in the case of the Ising model with a transverse field and transverse exchange, and in the case of macroscopic spins with uniaxial anisotropy. Our approach provides a transparent theory of level crossings, where the tunneling between time reversed configurations vanishes as a function of the external field. The crossings result from destructive quantum interferences between competing flipping channels. Our theory consistently predicts N crossings in chains of N Ising spins, 2S crossings in single spins of magnitude S, and yields explicit analytical formulae for the level crossings of open chains and macroscopic spins. Disorder can be easily implemented in this perturbative formalism. Leading disorder effects can be treated analytically for spin rings. At the smallest transverse field crossing, the suppression of tunneling is most robust with respect to disorder. We briefly discuss the implications of our findings for the use of realistic spin clusters on surfaces to store information.

Quantum tunneling in magnetic clusters has been intensively studied in the nineties as a special case of macroscopic quantum tunneling [1]. Quantum tunneling between two states with very different quantum numbers, e.g. S z = S and −S for large spin S, is in general a very high-order process since elementary terms such as a transverse field or exchange processes only change this quantum number by 1 or 2, and high-order perturbation calculations of the tunneling were limited to systems with a single tunneling channel [2][3][4]. Other approaches included a WKB approximation in the semiclassical limit [4,5], but the most successful approach proved to be an instanton and path integral formulation [6][7][8][9][10][11]. The predictions of these theories, for instance the difference between half-integer and integer spin [8,9], or the presence of oscillations of the tunneling as a function of a transverse field [10], have been beautifully confirmed by experiments on ferromagnetic molecules that measured Landau-Zener transition probabilities, which are sensitive to the tunneling between nearly degenerate levels [12]. The tunnel splitting was found to oscillate with the transverse field, and the position of the minima of the tunneling amplitudes were shown to alternate depending on whether the difference between the S z components of the initial and final spin states was even or odd. For a review, see Refs. [13,14].
More recently, it has become possible to create and control very small clusters of magnetic adatoms deposited on surfaces [15][16][17], where the exchange couplings between adatoms can be tailored by their positioning [18]. One promising idea is to use arrays of small antiferromagnetic (AF) chains or ladders of Ising spins as a means to store bit information in a very compact way [16]. Applying a large enough voltage pulse with the STM tip, tunneling between the AFM Ising ground states can be induced, providing a way to switch between bit states. The limitation of the associated memory comes from spontaneous tunneling, thermal or quantum, between the two AF ground states. For low enough temperature, the switching rate between the Ising AFM ground states saturates to the quantum tunneling rate, which decays exponentially with system size. Thus, one way to preserve the bit state longer is to increase the size of the cluster, at the cost of lower density of information. If however quantum tunneling depends strongly on an applied transverse field, with marked minima as in the case of molecular magnets, one could reduce the rate of quantum tunneling without increasing the cluster size.
The first indication that this could indeed be the case comes from a recent experiment that demonstrated that anisotropic chains of N spins-1/2 can have ground state crossings as a function of an applied magnetic field [19]. The crossings occur between the two lowest levels which form a quasi-degenerate subspace (which becomes ex-actly degenerate in the thermodynamic limit). At these crossings, the quantum tunneling between the two states is completely suppressed. For small spin chains, this can easily be demonstrated numerically [20]. A general theory of these level crossings has not yet been developed, however. One step in this direction has been achieved in Ref. [21]. This approach relies on a mapping of the spin model onto a fermionic chain using Jordan-Wigner transformation, and it is thus limited to open chains. A mean-field decoupling of this fermionic model maps it onto the Kitaev model of a one-dimensional p-wave superconductor [22]. This leads to the prediction of exactly N crossings, where N is the number of sites of the chain. In this framework, the crossings are naturally interpreted in terms of the oscillating, exponentially weak coupling between the two Majorana edge states on either end of the chain.
In this paper, we develop a perturbative approach that provides a unified and general framework for all these phenomena. It relies on a reformulation of degenerate perturbation theory to lowest non-trivial order in terms of a simple recurrence between flipping amplitudes. This approach leads to a general expression of the tunneling amplitude as a homogeneous polynomial in the amplitude of the transverse exchange processes and of the square of the transverse field. In this approach, the destructive quantum interferences that lead to level crossings appear as a natural consequence of the competition between tunnel processes with positive and negative amplitudes. In the case of a single exchange channel, level crossings will be present as soon as the amplitude of this process is positive (AFM) since it will then compete with the second order process due to the transverse field, which is intrinsically negative. Our approach leads to a number of analytical results in the limit of small exchange processes and transverse field (e.g. the exact solution for the open chain, or for a macroscopic spin), to very good asymptotic expressions for large closed rings and the effect of weak disorder, as well as to general qualitative conclusions (e.g. concerning the number of level crossings in a ring).
In higher dimension, our method still works, as we show on small rectangular, triangular, and cubic clusters, but it becomes rapidly very complex because of the large number of cluster shapes generated in the recurrence. The interest of the method is to show for the example of small clusters that the interference leads to a number of ground state crossings equal to the number of Ising spins, independently of the geometry, suggesting that this remains true for larger systems.
Our formalism turns out to be particularly convenient to study the effect of disorder (such as heterogeneities among different clusters, spatial g-factor variations, weak random fields, etc.) on the suppression of tunneling amplitudes. As we shall see, the suppression of tunneling in the lowest transverse fields is the least affected by disorder. This leads us to the conclusion that in order to achieve a maximally robust suppression of tunneling over a range of potentially fluctuating parameters of clusters, one should use the lowest transverse field for which a ground state crossing occurs in the disorder-free limit.
The paper is organized as follows: In Sec. I we introduce the spin models used throughout the rest of the paper. In Sec. II we present the iterative perturbation theory method applied to Ising models. In Sec. III, we apply the method to a 1D ring and chain, and in Sec. IV to small 2D and 3D clusters. In Sec. V we introduce disorder on the pertubative couplings of a ring and obtain the mean square displacement of the transverse field zeros, as well as the second moment of the tunneling to leading order in the disorder strength. In Sec. VI we apply our method to an anisotropic single spin model and obtain the tunneling amplitude, including an exact result for the crossing field values. In Sec. VII we discuss other systems where one can observe interference between different tunneling paths, and consequently crossings. Finally, in Sec. VIII we discuss our results and their experimental consequences.

I. MODELS
All the models we consider are spin models and can be written in the form where H 0 is a dominant diagonal term with a doublydegenerate ground state where the two states transform into each other by flipping all spins, and where V 1 and V 2 are perturbations that respectively flip one or two spins. λ is only an auxiliary parameter which we introduce to organize the perturbation expansion. It will be set to 1 later on. We thus require that the matrix elements of V 1,2 are much smaller than the norm of the terms in H 0 . We use the following notation for the two lowest energy states of H, Our main goal is to calculate in leading order in λ the energy splitting ∆E ≡ E + − E − . All the Hamiltonians we consider share the same symmetry, which physically corresponds to a spin reflection across the x-y plane, S z → −S z . We write the symmetry formally in a way that applies to all models, as where S z is the total spin projection along z, so that exp(iπS z ) rotates all spins by π around their z axis, and T is the time-reversal operator which inverts all spins. Let us call the two ground states of H 0 as |∅ , |Σ . The operator R flips all spins and transforms one ground state into the other, FIG. 1. Color plot of the energy difference between the two lowest energy levels of the model (7) on a spin chain of 8 spins, obtained with exact diagonalization. The full lines indicate the locations where the ground state is degenerate, while the dashed lines indicate our perturbative prediction (Eq. (61)), which works very well for small fluctuation parameters Jx and hx. Calculating the energy gap close to the phase transition hx = Jz, Jx = 0 perturbatively in the small parameters hx−Jz and Jx yields ∆ ≈ 2(hx − Jz − 2Jx) [23]. This gives the approximate location of the quantum critical line between the Ising-ordered and the paramagnetic phases as Jx ≈ (hx − Jz)/2, as shown in the plot. This is still accurate for large Jx since a similar calculation close to the phase transition at (hx = 2Jx, Jz = 0) gives the same gap [24]. For the full phase diagram, see Ref. [25]. There are no zeros when Jx < 0, and the ordered phase shrinks with increasingly negative Jx. Both reflect the fact that the transverse terms reinforce each other for Jx < 0, while they compete for Jx > 0.
Taking this into account together with the fact that R is anti-linear, we may write the (unnormalized) ground states as simultaneous eigenstates of H and R, as with two unknown phases φ ± .
In the following subsections we introduce the spin models used throughout the paper. We start with an Ising model with transverse terms, followed by other Ising models which are extensions of that model. One can solve the extended models straightforwardly, only requiring a transformation of the couplings to be able to use the solution of the first Ising model. Finally, we introduce an anisotropic single spin model, a conceptually simpler model than the Ising models, which was studied in the literature to analyze ground state degeneracies, having many similarities to a ferromagnetic Ising model.

A. Ising model with transverse field and exchange
We first consider the Ising model in a small transverse field and subject to weak transverse exchange interac-tions [21]: Here σ are the Pauli matrices and N is the number of Ising spins. We keep the lattice general for now, however, we do require that the lattice be classically unfrustrated, such that the ground state manifold of H 0 is only doubly degenerate. In the following sections we will carry out more detailed calculations by restricting ourselves to specific lattices and Ising interactions. We denote the Ising eigenstates and energies as where |m are the classical Ising configurations with the spins either up or down along the axis z. We label the states |m by the set m of spins that are flipped with respect to one of the Ising ground states. By definition, our reference ground state corresponds to the empty set m = ∅, and the other ground state to m = Σ, the set of all spins. The energy of these two states is ε ∅ (= ε Σ ). Before proceeding with introducing various extensions, let us comment on the ground state degeneracies on the spin chain with even N . When J x > 0 and for finite size, there are lines in the (J x , h x )-plane where the ground state is degenerate, as shown in Fig. 1. In our perturbative regime, these lines scale as h x ∼ J x |J z |. The lines continue for larger J x , h x , but the scaling becomes linear. In fact, when J x , h x J z , all lines except one approach the critical line J x = h x /2 of the J z = 0 classical model (an AFM Ising with longitudinal field), separating the AFM-ordered phase along x from the PM phase. The classical ground state on this critical line is highly degenerate, but the J z exchange coupling lifts the degeneracy and induces many ground state crossings close to the classical critical line. The line corresponding to the smallest field for given J x instead approaches the classical critical line J x = h x . On that line it costs no energy to flip the terminal spin of the chain which is anti-aligned to the external field in the AFM phase. This separate degeneracy line is present only for chains with even N .

B. Ising model with a staggered field
The first extension we consider is a staggering of the magnetic field on bipartite lattices with the same number of sites on either sublattice, but different transverse fields h x,A and h x,B acting on the two sublattices A and B, respectively. Thus, the modified perturbation is with unchanged V 2 . In Fig. 2  We will also extend the model by having transverse couplings along both transverse axes, On a 1D chain of spins with open boundary conditions, this model has an interesting limit in which it is exactly solvable. Taking J y = h x = 0, one can map the model with a Jordan-Wigner transformation [26] onto the Kitaev chain model [22], a free fermion model which one can solve exactly [21,27] to obtain field values in which not only the ground state is degenerate, but all eigenstates are: where n = 1, . . . , N/2 . If we take J x > 0, these zeros only appear if J z > 0, as one can see in Fig. 3. Interestingly, our perturbative approach yields exactly the same expression in the appropriate limit (Eqs. (61, 43)), with no higher order corrections. This is presumably a consequence of the fact that the zero energy Majorana fermions do not backscatter, so that our leading order approximation, which is equivalent to a forward scattering approximation for the fermions, becomes exact [28].

D. Single spin model
Instead of studying a cluster of N spins explicitly with all its internal couplings, one may consider the cluster as a big effective spin S and consider the tunneling from the up to the down state of this composite spin. The precise Hamiltonian for this equivalent big spin would in general be very complex. However, qualitative features can be expected to be captured by simple effective interactions which can be written as low powers of the total spin operators S x,y,z .
A transverse field on a ferromagnetic cluster of spins naturally translates into a transverse field acting on the big spin. Also, the ferromagnetic Ising configurations are the projections with largest total S z of the largest spin one can form with N (S = 1/2)-spins. Thus, the ground state of a Hamiltonian of the form −S 2 z corresponds to a ferromagnetic Ising ground state of the original cluster. Indeed, we shall see later that this model captures qualitatively the features of a ferromagnetic Ising model (J z < 0) as considered in Sec. I C. In contrast, it is not clear whether such an approximate mapping is meaningful for antiferromagnetic clusters.
Here we will reconsider quadratic single spin Hamiltonians equivalent to those considered earlier in the literature [3,5,10,12]. In particular, we take the most general quadratic single spin model with anisotropy and a field transverse to the easy axis. Upon choosing axes that bring the quadratic part to a diagonal form, we are left with three quadratic couplings J x,y,z . We are free to choose the easy, medium and hard axes to be, respectively, along the z-, y-and x-axes (i.e., J z ≤ J y ≤ J x ). We are also free to set one of these couplings to zero, since the ground state splitting will be independent of a constant term J S· S that one can add to the Hamiltonian. We therefore consider the Hamiltonian where we chose J y = 0, implying that J z ≤ 0 and J x ≥ 0.
When performing perturbative calculations, we assume that This model has 2S = N values of the magnetic field where the ground state is doubly degenerate, as demonstrated by A. Garg [10] using path integral and instanton calculations justified at large S, but without restricting to the perturbative regime (13). In the considered large S limit, the crossings can be interpreted as the negative interference of two tunneling instantons having different Berry phases. We shall derive a very similar result within perturbation theory, but without relying on the size of the spin (or N ) and without taking a saddle point approximation. We will show in Sec. VI that here again the underlying mechanism behind the zeros is the competition of multiple tunneling paths with oscillating signs. Finally, using a different approach we will derive non-perturbatively the location of the 2S equally spaced transverse field zeros without relying on any approximation.

II. PERTURBATION THEORY FOR COLLECTIVE TUNNELING
In this section we will present the perturbative method applied to Ising models defined in Sec. I. The method is derived in more detail in appendix A.
Standard techniques to carry out degenerate perturbation theory at arbitrarily high orders construct a perturbative expansion for an effective Hamiltonian H eff which only acts on the Hilbert space spanned by unperturbed states and yields the exact splitting of the ground state manifold due to perturbations. Let us call the set of unperturbed ground state configurations g = {∅, Σ}. Defining P as the projector onto the subspace g, one constructs H eff which projects out all the excited states, i.e., P H eff = H eff P = H eff . However, H eff is not fully specified by these requirements. A first full series expansion for a possible choice of H eff was obtained in [29] for a general Hamiltonian with a discrete spectrum. It leads to a generalized eigenvalue equation that must be solved for the split ground state energies. A variation of [29] was later given in [30], where the eigenvalue equations are simpler, but the operator H eff will in general turn out to be non-Hermitian. Here we use the latter approach. The eigenvalue equations of this effective Hamiltonian are where |Ψ ± are the lowest energy eigenstates of H. Owing to the anti-linear symmetry R (Eq. (4)), which lets us write |Ψ ± as in Eq. (6), we deduce that the energy splitting is given by in terms of the off-diagonal matrix element of the effective Hamiltonian. Our perturbative method allows us calculate the matrix element ∅| H eff |Σ to leading order in λ. Since V 1 and V 2 respectively flip one or two spins when acting on |m , and given that V = λV 1 + λ 2 V 2 , the power of λ of a tunneling path between |∅ and |Σ corresponds to the number of spin flips that occurred. Since the shortest paths have exactly N spin flips (each spin flips once and only once), it follows that λ N is the lowest order that will occur, and thus where we have defined t N as the leading order term.
In the limit λ → 0, the phases in Eq. (6) must vanish, φ ± → 0. This follows from the fact that in this limit the ground states approach the two linear combinations One thus finds As shown in Section A, the tunneling amplitude is given by where the l i=1,...n ∈ {1, 2} obey l 1 + l 2 + · · · + l n = N, and where ∆ε m ≡ ε m − ε ∅ is the unperturbed excitation energy of the state |m . Let us denote by |m| the cardinality of the set m, that is, the number of spins that are flipped relative to the ground state ∅, so that |Σ| = N . The tunneling t N can be calculated recursively. To this end we introduce c m as an intermediate tunneling coefficient, analogous to t N in Eq. (19), where we impose the condition l 1 +l 2 +. . . l n = |m|, on the l i , which restricts the sum to the terms that contribute to leading order (λ |m| ) to the tunneling between |∅ and |m . Finally, by inserting the identity as m |m m | before the last two factors and expanding V ln S |m , one obtains a recursion relation that connects c m to coefficients of smaller clusters m ⊂ m, yielding the cluster recursion relations where we sum over flipped spins i ∈ m, or pairs of flipped spins {i, j} ∈ m. An analogous expression follows for t N , The recursions (22,23) hold for general perturbations V 1,2 which flip single spin or pairs of spins, respectively. It is straightforward to generalize the recursions (22) to include higher order terms of the form λ k V k , which flip k spins. The iterative procedure can be further simplified due to the independence of disconnected clusters. Let us consider a cluster m that is composed of several mutually disconnected clusters m i of spins flipped relative to ∅. We call a set of clusters disconnected if the excitation energy of the set is the sum of the excitation energies of the individual clusters, i.e. ∆ε m = ∆ε m1 + ∆ε m2 + . . . ∆ε mn . In Sec. B, we show that for such separable sets, the intermediate tunneling coefficient is the product of tunneling coefficients for their disconnected components, To calculate t N for a given cluster one proceeds as follows: One identifies all inequivalent connected subclusters m of the considered spin cluster, associates a coefficient c m to each of them, and uses Eqs. (22,24) to calculate the coefficients recursively for increasing cluster sizes |m|. We say that two clusters are equivalent if their coefficients c m are the same, which is for instance the case if the two clusters are symmetry related. If the unperturbed Hamiltonian is the Ising model with nearest neighbor couplings, the excitation energy ∆ε m of a cluster is 2|J z | times the number of bonds that connect flipped spins (∈ m) to unflipped spins (∈ Σ\m). In other words, the excitation energy is proportional to the total length of all domain walls between m and its complement.
In general, the recursions (22,23) can be rather complicated to solve, especially if one has to consider a large number of inequivalent clusters. However, we will see in the following subsections that they take a simplified form when applied to some of our models, and even result in recursion relations with closed form solutions in some cases. Applying the method to the first Ising model (7), our next step, will be very instructive.

A. Tunneling in Ising models with transverse field and exchange
Let us consider the model (7) and let us apply the recursion relations (22,23) to it. The first thing to note is that the matrix elements in the recursion relations simplify greatly. They are given by and if i, j are nearest neighbors, and vanish otherwise. Then, the recursion relations become In this model, two clusters are equivalent if they are identical including their environment up to their first neighbors. As an example of this, consider a cluster of flipped spins in the bulk of a lattice with open boundary conditions. Any translation by a lattice unit whereby the cluster does not touch the boundaries results in an equivalent cluster. When a boundary is reached instead, the first neighbors of the cluster change and we will find a different cluster coefficient.
Since t N is the sum over all tunnel paths that flip every spin exactly once, it is clear that the resulting expression is a polynomial in h x and J x , each term being proportional to a product J k x h N −2k x with 0 ≤ k ≤ N/2 . The general form for the tunneling coefficient of a cluster of N spins, regardless of the lattice, will thus take the form with some lattice-dependent real coefficients a k , and (potentially complex) roots α n . Due to the minus sign in the projector S (Eq. (20)), the sign of a k reflects the number of flipping terms V 1,2 that are applied on the corresponding tunneling paths. The sign thus alternates with k: where we have also taken into account the negative sign of the matrix element of V 1 in Eq. (25). If J x < 0, all monomials contribute with the same sign, and thus |t N | grows monotonously with the field h x , with no zero crossing. However, if J x > 0, there is a negative interference between paths with different numbers of perturbative steps, and t N may oscillate as a function of the field. In this case, we can have ground state crossings. If the polynomial of Eq. (30) has positive real roots α n , the crossings occur at the fields Thus there may be up to N values of the transverse field (or N/2 , if one restricts to h x ≥ 0) where t N = 0, depending on the number of real α n . (Note that real α n are necessarily positive, since negative α n would imply zeros for J x > 0, which is excluded).
Later on we will solve the recursion relations for specific spin clusters and lattices, where we do find that in all cases considered the α n are real and positive. For now, let us assume that indeed all α n > 0, so that we have either N degeneracy points or none (except for the trivial one h x = 0 for odd N ), depending on the sign of J x .

B. Tunneling in Ising models with a staggered field
We consider the model presented in Section I B where the field of the model (7) is staggered. One could calculate the matrix elements and write down the recursion relations for this model, but this turns out to be unnecessary. Let us first guess the polynomial form of t N for this model. The tunneling paths due to transverse exchange only contribute with a factor J N/2 x . For tunneling paths involving spins flipped by transverse fields, those must come in equal numbers on the two sublattices, and thus t N must take the form where the coefficients a k must be the same as those of the polynomial for a homogeneous field h x , cf. Eq. (29). Thus one simply should substitute h 2 x → h x,A h x,B in that equation. It follows that ground state degeneracy occurs whenever The condition h x,A h x,B J x > 0 is a general prerequisite for such degeneracies. If all α n are positive, we recover the behavior observed numerically in Fig. 2 C. Tunneling in Ising models with general transverse couplings In the case of the more general Ising models (10), we found the behavior of the ground state crossings to depend on the Ising ground state (Fig. 3). In particular, we distinguish whether the dominant Ising interactions are ferromagnetic (J z < 0) or antiferromagnetic (J z > 0), respectively.

Ising ferromagnets
Since the ground state is ferromagnetic, and since we only flip each spin once, the matrix element m \ {i, j}|V 2 |m appears in the recursion relations only in the form for two neighboring spins, while m \ {i}|V 1 |m only appears as This implies that the resulting recursion and the tunneling amplitude t N will be the same as Eqs. (27,29), up to the substitution and it suffices to solve the model (7). Performing the substitution in Eq. (32), one finds that the tunneling only vanishes if either h y = 0 and J x > J y , in which case there are zeros at the fields or if h x = 0 and J y > J x , at fields In other words, the transverse field has to be applied in the spin direction which corresponds to the stronger antiferromagnetic (or weaker ferromagnetic) exchange. This is analogous to the result found in the single spin case by Garg [10], which we will review in Section VI below.

Ising antiferromagnets
For antiferromagnetic Ising models case, let us consider a bipartite lattice, so that the ground states have opposite spins on each sublattice. V 2 only acts in the form while V 1 flips spins from down to up on one sublattice and from up to down on the other sublattice, We take each sublattice to have the same number of spins.
Since the lattice is bipartite and V 2 flips one spin from each sublattice, there must be an equal number of single flips due to V 1 on either sublattice. This implies that we can simply substitute in Eqs. (27,29) to obtain the result for the generic, bipartite antiferromagnetic Ising models. Interestingly, the direction of the homogeneous transverse field in the x − y plane is irrelevant. That is, zeros of the tunneling amplitude occur in circles in the transverse field plane, provided that the transverse exchange is predominantly antiferromagnetic (J x + J y > 0). The tunneling vanishes for transverse fields of magnitude regardless of its angle in the x − y plane. While there is no angle dependence to leading order, the radial symmetry is broken by higher order corrections as confirmed numerically in Fig. 3. Note that the tunneling matrix element t N is in general complex, and thus the condition t N = 0 determines a manifold of codimension 2 in the parameter space of transverse couplings. Thus, by fixing the exchange couplings J x , J y and looking for zeros in the transverse field plane, one will generically only find isolated points, as it happens in the case of ferromagnetic clusters. A qualitatively different situation arises in antiferromagnetic clusters because there, owing to the substitution (42), the tunneling amplitude is always real, such that zeros organize in a manifold of codimension 1, i.e., closed lines in the transverse field plane.

III. TUNNELING IN 1D SYSTEMS
In this section we apply our method to the 1D model (7). The recursion relations for 1D systems are rather simple because any connected cluster of flipped spins is uniquely defined by its length and position. We first consider a ring of N spins, where the exact solution of the recursion allows us to extract explicit asymptotic expressions for large N . Then we consider an open chain, where we even obtain a closed analytical expression for t N for any N .

A. Closed chain: a ring of spins
For a ring of spins, all connected clusters of a given length are equivalent. We denote by c l the intermediate tunneling coefficient associated with a cluster of length l. The application of Eq. (27) is straightforward: By unflipping a single spin from the cluster, we obtain one of l possible states with l−1 flipped spins. If the unflipped spin is at the edge of the cluster, its contribution to c l is h x c l−1 /(4|J z |), since the excitation energy of the cluster is 4|J z | due to the two domain walls at its ends. If the unflipped spin is in the bulk then the state consists of two clusters. We thus use Eq. (24) to write the contribution of that state as h x c n c l−n−1 for some 0 < n < l. Defining c 0 ≡ 1 we can combine the edge and the bulk terms, and using an analogous reasoning for the term related to J x we finally have To obtain the polynomial for the full tunneling coefficient t N , a similar recursion can be used. It slightly differs from the above due to the periodic boundary conditions. Unflipping a pair or a single spin in the ring, we are left with a single cluster of length N − 2 or N − 1. The unflipped spin(s) can be at N positions, so that we find In general the location of the field-zeros h (n) x depends on N . Interestingly, it turns out that the pair of largest zeros, ±h , is common to chains of any size and takes the value h To show that this is indeed so, we start from Eq. (45). The condition to have t l+1 = 0 requires c l = (J x /h x )c l−1 . If this is to hold for all l and given that c 0 = 1, we must have c l = (J x /h x ) l . Using this in Eq. (44) and simplifying we find that this relation is indeed satisfied if h x = ±2 J x |J z |. For this value of h x the tunneling t N thus vanishes for any N . We will retrieve this result from a direct calculation of t N below.
In order to calculate t N for any value of h x , we define the generating functions where z is a complex variable. We then multiply Eqs. (44, 45) by z l , and sum them from l = 2 to l = ∞. Solving for C(z) [31] and eventually for T (z), we obtain the closed expression where The above formula for T (z) represents the power series (47) with its domain of convergence at small enough z, but analytically continues it beyond. We can now calculate t N by contour integration of T (z)/z N +1 around its pole at z = 0. We have the exact expression We now deform the contour, pushing it to infinity, but avoiding the branch cuts ending at z ± . This is best done using keyhole contours (Fig. 4). The precise contour used depends on whether z ± are real or a pair of complex conjugate numbers, which we discuss separately. We restrict ourselves to J x > 0, since only in that case t N exhibits interesting oscillations. a. h 2 x > 4J x |J z | -In this case, both z ± are real and either both positive or both negative, depending on the sign of h x . Without loss of generality, we take h x > 0. The appropriate contour is shown in Fig. 4 on the left. Pushing the radius of the large circle to infinity and letting the radius of the small circles around z ± shrink to zero, the integrals 1 and 4 cancel, while the integrals 2 and 3 between the branch points add up to Note that upon changing h x within the domain h 2 x > 4J x |J z |, the integrand remains positive, and thus t N never becomes zero.
b. h 2 x = 4J x |J z | -At the border of the above domain one has z + = z − . From Eq. (50) one sees that the tunneling becomes zero at this point, independently of N , as we have already found previously. Since there are no zeros at higher fields, this corresponds to the largest field zero.
c. h 2 x < 4J x |J z | -Here the branch points become a pair of complex conjugate numbers. We have where We consider the contour shown in Fig. 4 on the right, where the branch points z ± are avoided with keyholes oriented radially along the lines z = re ±iθ0 . Shrinking the small circles to zero, and expanding the large circle to infinity, the expression for t N simplifies to the contribution of two line integrals along the radial branch cuts, resulting in the exact expression: (53) At large N , we can make progress by replacing re iθ0 − r 0 e −iθ0 by its value at r = r 0 (which is valid as long as θ 0 1/N ). The remaining integral can be written in terms of Gamma functions. To leading order at large N one obtains (54) Note that this expression has, however, N + 2 zeros as a function of h x : N zeros arise from the vanishing of the high frequency sine at fields given by with n = 1, . . . , N/2 . Two further zeros are due to the vanishing of sin θ 0 . Those reproduce correctly the pair of largest field zeros, h x = ±2 J x |J z |, which we have already identified above. The two zeros h (n) x with n = N/2 are instead artefacts that are introduced by approximating the numerator in the integrand with its value at r = r 0 . This approximation is not controlled in that field regime since there one has θ 0 = O(1/N ). These two zeros thus have to be discarded, and we are left with N zeros, as it should be.
In Fig. 5 we compare the asymptotic Eq. (54) with the exact polynomial for t N obtained from explicitly solving Eqs. (44, 45). The agreement is very good even for moderate N , and it further improves with system size.

IV. 2D AND 3D CLUSTERS
While in 1D chains connected clusters only come in one shape (a connected stretch of spins), the length and  Table I for which we have to compute the intermediate tunneling coefficients. The 3 × 2 cluster differs from a ring of 6 spins only by one additional bond. By showing that the tunneling polynomial still has N zeros, we thus demonstrate the robustness of the number of zeros to certain perturbations. Upon gradually turning on the bond that transforms the ring into the 3 × 2 cluster, the zeros move towards higher fields, as one expects. We further derive the tunneling polynomials t ∆ for an equilateral triangle made from 6 spins, cf. Table II, and t 2×2×2 for a cube of 8 spins, cf. Table III. Solving the resulting recursion relations given in the Tables I, II and III, we obtain the following polynomials, where we define the variable which have the expected number of zeros, N = 6 or 8, as one can see in Fig. 6.

V. WEAK DISORDER
In this section we consider a 1D ring with weak disorder in the transverse fields and in the exchange in the form of where J x,i connects spins i and i + 1 and where J x > 0, such that there are zeros in the absence of disorder. We denote the disordered tunneling ast N , reserving t N for the tunneling in the disorder free limit. One can in principle determine the recursion relations of such a system and thus study disorder using explicit polynomialst N . However, we shall rather calculate the relevant average quantities to lowest order in an expansion in δh x,i and δJ x,i . A quantity of particular interest is the typical finite tunneling induced by random fluctuations of the couplings when the average external field is held at one of the tunneling zeros h (n) x . Let us denote the disorder induced tunneling at a zero as T n ≡t N hx=h (n) x . (66) We are also interested in how much a ground state crossing shifts due to the presence of randomness: whereh (n) x and h (n) x are respectively the n th zeros of the polynomialst N and t N . We can calculate the second moment of these random variables by considering a homogeneous ring where all spins and all bonds are equivalent and where in the absence of disorder all couplings are equal. However, explicit calculations of the second moment of T n and K n for open spin chains showed qualitatively very similar behavior as we find below for rings. For a ring, to first order in the perturbations, one has This follows since perturbations on different sites are equivalent, and thus all partial derivatives are equal: and likewise for the exchange. For identically and independently distributed local disorder, the disorderinduced variance of the tunneling evaluated at a transverse field zero h (n) x thus results as To leading order, the response of the tunneling amplitude T n to disorder is linear. The shift of the transverse field, K n , necessary to compensate for this disorder-induced tunneling is then given by the relation From Eqs. (70, 71), we see that in order to calculate K n we only need the ratio of the derivatives . This ratio is easily obtained from the factorized polynomial form of t N (Eq. (30)) as From this we deduce the average mean square of the drift in the zeros to leading order as: where quantify the relative fluctuations of the couplings. In many cases these are the most appropriate measures of the disorder strength. We see that the zeros corresponding to larger fields are more strongly affected by disorder. This certainly holds as long as the shifts are smaller than the typical spacing O(1/N ) between zeros.
One can obtain the standard deviation of the tunneling amplitude from Eqs. (71, 74), calculating ∂ t N ∂h x from the previously obtained expressions. For a ring (Eq. (60)) we get Again, we see that T n ∼ h (n) x for small transverse field, while it grows quickly as the largest zero h

VI. TUNNELING IN SINGLE SPIN MODELS
Here we consider the single spin models as presented in Sec. I D. Similarly to an Ising ferromagnet, the unperturbed (λ = 0) ground states correspond to the two states with S z = ±S. Applying our method to this Hamiltonian is rather straightforward: We simply calculate the tunneling matrix element from −S to the +S ground state to lowest order in λ, using intermediate tunneling coefficients c m , where m refers to the spin projection onto the z-axis, with eigenstates defined by If we use the Hamiltonian in the form of Eq. (12), the recursion for c m will involve both c m−1 and c m−2 , but unlike in the problem treated in Sec. II A, the matrix elements and the denominators involved in the recursion depend themselves non-trivially on m. The resulting recursion is hard to solve analytically. We can, however, simplify the recursion greatly by first performing a rotation in the x − z plane: where we choose α to satisfy tan 2 α = λ 2 J x /|J z |, such as to kill the matrix elements k|H|k − 2 between S zeigenstates, S z |k = k |k . This yields the Hamiltonian in the rotated basis (79) where we consider dominant ferromagnetic Ising coupling J z < 0. In the perturbative regime, the rotation angle α is small and the ground states of H still have a large overlap with the two S z eigenstates |±S . We now deal with the problem of calculating the tunneling between these two states up to order λ 2S . Since the matrix form of H in the S z basis is tridiagonal, only the off-diagonal terms proportional to λ contribute to t 2S . Thus we write the Hamiltonian in the form where The matrix elements of these operators are where we used k|S x |k − 1 = −i k|S y |k − 1 , and Now that the Hamiltonian is tridiagonal, there is a single tunneling path between |±S of order λ 2S . Its contribution to t 2S is just the product of all off-diagonal matrix elements k|V 1 |k − 1 divided by (minus) the energies of all intermediate excited states. The total tunneling amplitude is .
Note that this already takes the form of a factorized polynomial, with zeros given by the equation where n = −S + 1, −S + 2, . . . , S. Recall that we took the x−axis to be the hard axis (J x > 0). We find 2S equally spaced ground state crossings for h y = 0 and transverse fields This is the same qualitative behavior as in a ferromagnetic Ising cluster as seen in Sec. I C, cf. Eqs. (38,39). This result coincides in lowest order with the zero positions as calculated by A. Garg [10] and as demonstrated experimentally by Wernsdorfer [12]. However, in the rotated frame it becomes quite simple to obtain the exact degeneracies of this model, as we show in the next subsection.

Exact degeneracies of single spin model
If we are only interested in the location of the zeros, the rotated Hamiltonian can be used beyond perturbation theory to determine the position of the zeros exactly. Indeed, a zero of the tunneling matrix element S|H | − S occurs whenever one of the off-diagonal entries of H becomes zero. At that point, the Hamiltonian splits into two uncoupled blocks for S z ≥ k and S z ≤ k − 1, implying that the up-state is strictly decoupled from the down-state to all orders. This entails an exact double degeneracy of the ground state (which map onto each other upon rotation by π around the x-axis). Now, the off-diagonal matrix element is given by Since J x > 0, we can again have a ground state degeneracy only for a transverse field along x. The critical fields are determined by the exact condition which agrees with the perturbative result of Eq. (89) to lowest order in J x , and reproduces the non-perturbative path integral results by Garg [10]. Here we have shown that this yields the location of the zeros exactly, independently of the size S of the spin. The result that the transverse field has to be applied along the hard axis is fully consistent with what we found for FM spin clusters in Sec. I C. The main difference between the exact cluster calculation and the single spin model consists, however, in the precise location of the zeros. For the single spin model, the zeros are equally spaced, while for clusters they are spaced more and more densely the larger the transverse field, as one can see, e.g. in Fig. 5, or read off from the analytical result in Eq. (55).

VII. OTHER SYSTEMS WITH COMPETING TUNNELING CHANNELS
The mechanism we have studied here, namely the interference of parallel multi-step tunneling channels between an initial and a final state is very general in nature and appears in various physical contexts.
A famous example is the case of resonant singleparticle tunneling via several intermediate sites, a problem introduced by Nguyen, Spivak and Shklovskii [32,33], with comprehensive reviews given in Refs. [34,35].
For free particles (non-interacting fermions), different paths from an initial to a final site contribute with an amplitude whose sign alternates with the number of intermediate sites whose energy is above the chemical potential. This leads to negative interference between alternative paths. A magnetic field introduces additional Aharonov-Bohm phases and decreases the likelihood of full negative interference, resulting in increased transmission, that is, negative magnetoresistance [36,37]. Recently, it was found that the equivalent question for hard core interacting bosons leads to a similar interference problem, where, however, at energies close to the chemical potential all path amplitudes contribute with the same sign, leading to maximally constructive interference. [38,39] This situation resembles that of a ferromagnetic cluster with transverse field applied in the direction in which the transverse exchange is more ferromagnetic (i.e., the softer axis). In contrast to the magnetic clusters, however, in these hopping problems it is very hard or even impossible to tune a parameter (e.g., the magnetic field or the chemical potential) to suppress the tunneling completely.
Competing tunneling terms also arise in more general magnetic clusters composed of electronic and nuclear spins, a situation that frequently occurs in rare earth compounds. The magnetic ions are coupled to their nuclear spins, while the electronic spins couple to each other via dipolar and/or exchange couplings. Clusters of such ions often have doubly degenerate ground or excited states, which are only split by higher order tunneling processes that involve the interference of transverse fields, exchange/dipolar interactions and hyperfine couplings, that generically contribute with competing signs. Tuning the transverse field often allows to induce zeros in the corresponding collective tunneling. Similarly, the tunneling of the spin associated with a crystal field doublet of a magnetic ion can under certain circumstances be suppressed by a transverse field applied at specific angles, if different channels involving the magnetic field and transverse crystal field terms compete.
Ground state crossings have also been reported in SU(2) invariant, gapped frustrated spin chains [40,41]. In that case, the crossings are related to the interaction between the edge states of the chain. This is reminiscent of the explanation of the level crossings in the model of Eq. (7) in terms of Majorana edge states [21], and it is natural to ask whether these crossings can also be seen as a consequence of destructive interferences between different channels. For that purpose, let us consider the level crossings in the bilinear-biquadratic spin-1 chain [41], If one adds a strong uniaxial anisotropy along z, one may work with respect to an AF ground state, and the transverse terms with ∆S z = 2 and 4 have competing signs if J b > 0, presumably leading to level crossings similar to those of the isotropic case. It would be interesting to see if a more direct connection can be established by studying the effective coupling between the edge states starting from the AKLT model J b = J 1 /3 for which the edge states are fully decoupled in the ground state [42]. This goes beyond the scope of the present paper, however.

VIII. SUMMARY AND OUTLOOK
High order degenerate perturbation theory allows us to understand transverse field zeros in terms of negatively interfering tunneling paths, which in turn is tied to the presence of competing quantum fluctuations in the Hamiltonian. Our method nicely applies to 1D systems, where the tunneling can be obtained for any system size exactly, in contrast with 2D and 3D clusters where the number of different connected clusters grows exponentially with system size. Overall, the results support the existence of N zeros in some region of the parameter space independently of the geometry.
The original model (Eq. (7)) can be further extended by staggering the field or by adding exchange couplings along the y-axis while keeping the crossings. As we saw in Sec. I C, systems with FM ground states exhibit zeros only when the field is applied along the "hard axis" (the one with the strongest antiferromagnetic or the weakest ferromagnetic coupling). In contrast, AFM clusters on a bipartite lattice exhibit suppressed tunneling on approxi-mate circles in the transverse field plane. This may make AFM cluster ground states more attractive since the tunneling suppression is more resistant to fluctuations in the orientation of the applied field. This ability to control and suppress the quantum fluctuations in small magnetic clusters or single molecule magnets is indeed considered an important goal [43].
Introducing disorder in the exchange couplings and in the transverse fields, the crossings change position but do not disappear. The latter only happens when in ferromagnets J x − J y , or in antiferromagnets J x + J y , starts to change sign and turn negative. The relation between a set of J x , either randomly generated or carefully chosen, and the resulting number of crossings remains to be studied more deeply.
A certain amount of disorder in the exchange is always to be expected from static sources such as lattice imperfections, strain, or dynamically due to slow phonons. Spatial inhomogeneities can also induce g-factor variations that lead to an effective disorder in the transverse field. In an ensemble of weakly disordered clusters it is thus impossible to suppress the tunneling simultaneously in all clusters, and even in a single cluster temporal fluctuations of the parameters will destroy the perfect negative interference of competing tunneling channels. The best strategy to suppress the tunneling as much as possible consists then in tuning the transverse field to the first (smallest) zero, h (n=1) x , corresponding to the average exchange coupling in the system. The disorder-induced fluctuations away from vanishing tunneling turn out to be smallest under those conditions. This is closely related with the fact that the location of this smallest transverse field zero moves the least as the parameters of the Hamiltonian are slightly perturbed. Hence this zero seems to be the most interesting one for most applications.
Our recursive calculation of collective tunneling amplitudes generalizes nicely to simpler single spin models, and the ground state crossings in this model can be interpreted with the same tunneling interference argument. Given the 2S crossings of a single spin, one may expect that an appropriately chosen spin-S model on a lattice of N spins will exhibit 2SN crossings.
It is an interesting question to ask what happens to the zeros as one leaves the perturbative regime. In the ferromagnetic single spin model we can trace them easily, since we can obtain them exactly. If the hard axis is along the x-axis and one tunes J y up to and beyond J z for example, the number of 2S zeros remains intact, even though the easy axis has undergone a flop from the zto the y-axis. In lattice models, the zero lines in the h x -J x plane do not seem to disappear either. They even may cross quantum phase transition lines, as long as they enter a new phase with a degenerate ground state or a gapless phase. The study of the related phenomena and implications is left for future work.
Appendix A: Deduction of method Following Bloch's recipe [30] we consider the effective Hamiltonian H eff projected onto the unperturbed ground state subspace g = {∅, Σ} by the projection operator P . It takes the form where n specifies the order in V of the term. For a given n, we sum over all (n − 1)−tuples of k i = 0, 1, . . . that obey for all s = 1, . . . , n − 2, and k 1 + k 2 + · · · + k n−1 = n − 1.
The operator S k is defined as where ∆ε m = ε m − ε ∅ . The eigenvalue equations read where |Ψ ± are the lowest energy eigenstates of H. Due to the symmetry R (Eq. (4)) we may write the eigenstate projections to leading order as P |Ψ ± = |∅ ± |Σ + O(λ).
The relevant matrix element to calculate is Upon expanding V = λV 1 + λ 2 V 2 in Eq. (A1) and substituting in t, we group terms according to their powers of λ. Taking n 1 (n 2 ) to be the number of V 1 (V 2 ) operators present in a term of orders n and u, we have n 1 + n 2 = n, n 1 + 2n 2 = u, from which we find what orders of n contribute to u by taking the limiting cases of n 1 = mod(u, 2) and n 2 = 0. Finally, we only need to sum over the permutations of V 1 and V 2 that respect the order u. Applying this to t, we have where the l i = 1, 2 obey l 1 + l 2 + . . . l n = N. It helps to look at the calculation of the matrix element in Eq. (A9) sequentially; that is, starting with the extreme left operator, we apply each operator to the states on its left. While V 1 and V 2 always transform the states m| they act on, S k mainly acts as a projector onto a subspace of states, either the ground state g (if k = 0), or the excited states (if k > 0). Now, the treatment of V 2 as a second order perturbation is crucial for our method and physically justified by the fact that the basic action of V 2 on |m is to flip pairs of neighboring spins while V 1 flips single spins. In the term of order λ u , ∅| is acted on with enough V 1 's and V 2 's to at most flip u spins. Since we need a minimum of N spin flips to transform ∅ into Σ, it follows that the lowest order is λ N . After applying a V l to the states on its left, the resulting states must have l more spin(s) flipped than before for such terms to yield a non-zero contribution to order λ N . In particular, this means a projection onto g by S 0 = −P would only give terms that will eventually have zero contribution. This imposes k i > 0 in the leading term O(λ N ). However, the constraint (A3) only allows for one single choice of the k i , namely, k i = 1 for all i. Using this information in Eq. (A9), and writing t = t N λ N + O(λ N +1 ), we now have and we remind the reader that |m| is the number of spins of |m that are flipped relative to |∅ . Summing over l n and redefining n → n − 1, we have where l 1 + . . . l n = |m| − 1, l 1 + . . . l n = |m| − 2.
By expanding V 1,2 S |m in Eq. (A13), we shall see that we recover the cluster coefficients of smaller clusters. Consider first V 1 S |m . We have that which follows from the definition of S. Then, applying V 1 to the state |m , we get a sum over states m which differ by one spin flip from m. However, only m clusters with |m | = |m| − 1 yield a non-zero contribution to c m . Thus, With an analogous argument applied to V 2 S |m we obtain the recursion relations (23,22) in the main text. The expression for t N follows from an analogous derivation, the only difference being that there is no insertion of the operator S at the last step, as one can note from comparing Eq. (A10) and Eq. (A11). This eliminates the corresponding energy denominator.

Appendix B: Cluster independence
Consider a cluster C which is composed of two (dis)connected clusters of flipped spins A and B, by which we mean that the excitation energy of cluster C is the sum of independent excitation energies, Here we want to prove the relation among intermediate tunneling coefficients: which one should expect to hold because to leading order we can simply reduce the Hamiltonian to the parts acting on either A or B and drop all other terms, so that the flipping of A and B are independent processes. This assertion is conveniently proved by induction on the size of the cluster C. We suppose that we have proved it for a small cluster of size up to |C| − 1. For size 0 and 1 the assertion is trivial. The general recursion formula (22) shows that where we sum over the m obeying |m | = |C| − 1 or |m | = |C| − 2. Let us now write C = A ∪ B, and m = A ∪B or m = A∪B , depending on where V acts. Thus: where in the second line we used the induction hypothesis for smaller clusters, which implies that c A ∪B = c A c B . Now we use the relation (B3) in the form and an analogous expression for c B . Together with (B1), this proves the relation (B2).