Exact resummation of the Holstein-Primakoff expansion and differential equation approach to operator square-roots

Operator square-roots are ubiquitous in theoretical physics. They appear, for example, in the Holstein-Primakoff representation of spin operators and in the Klein-Gordon equation. Often the use of a perturbative expansion is the only recourse when dealing with them. In this work we show that under certain conditions differential equations can be derived which can be used to find perturbatively inaccessible approximations to operator square-roots. Specifically, for the number operator $\hat n=\hat a^\dag a$ we show that the square-root $\sqrt{\hat n}$ near $\hat n=0$ can be approximated by a polynomial in $\hat n$. This result is unexpected because a Taylor expansion fails. A polynomial expression in $\hat n$ is possible because $\hat n$ is an operator, and its constituents $a$ and $a^\dag$ have a non-trivial commutator $[a,a^\dag]=1$ and do not behave as scalars. We apply our approach to the zero mass Klein-Gordon Hamiltonian in a constant magnetic field, and as a main application, the Holstein-Primakoff representation of spin operators, where we are able to find new expressions that are polynomial in bosonic operators. We prove that these new expressions exactly reproduce spin operators. Our expressions are manifestly Hermitian, which offer an advantage over other methods, such as the Dyson-Maleev representation.


II. INTRODUCTION
Square-roots of operators appear in a large number of contexts in theoretical physics, and also play an important role in operator theory. In some cases, it is practical to calculate the operator square-root (OSR) using explicit formulas or by diagonalizing the operator. Often, however, there is only a very limited set of analytical tools to treat them, typically in the form of perturbative expansions. This is not because OSRs represent a niche problem. Indeed one of the earliest appearances was near the beginning of quantum mechanics in the square-root of the Klein-Gordon equation [1][2][3][4][5][6][7]. Even for such an old problem it may prove useful to have a larger analytical toolbox. Another prominent example of OSRs occurs in the Holstein-Primakoff spin representation [8,9], which is the usual starting point for spin-wave theory calculations. A third important OSR shows up in the context of quantum information in the form of the fidelity function [10][11][12][13][14][15][16][17][18][19][20][21][22][23] and the Bures metric [12,17,20,[23][24][25], both used to quantify the closeness of two quantum states. The purpose of the current paper is, however, not to review all examples of OSRs, but to introduce a non-perturbative approximation of OSRs.
Our method is inspired by several flow equation approaches to many-body problems. For instance, the Wegner flow equation approach [26], which was applied to various problems [26][27][28][29][30][31][32][33][34][35][36][37][38][39], allows for a non-perturbative diagonalization of a Hamiltonian using flow equations for its couplings. In this approach, the problem of diagonalization is recast in terms of differential equations. Similar methods have recently been used by some of us to find effective Floquet Hamiltonians [40], and various approximations to the time evolution operator [41]. Differential equation approaches have also been used in the method of unitary integration for the Liouville-Bloch equation [42] and Lindblad equation [43]. We aim to use a similar approach to approximate an operator square-root.
The application that may be of most current interest is the Holstein-Primakoff (HP) OSR. The HP representation is typically used in the context of spin (local moment) models to represent deviations around a welldefined spin order in terms of a single species of boson per lattice site. It allows for a perturbative expansion in the number operator of such bosons, and ultimately leads to linear [44] and non-linear [45] spinwave descriptions of quantum magnets [9]. However, for many systems of in-terest a ground-state spin ordering may be unknown or fail to exist, such as in frustrated systems [46], spin liquids [47][48][49][50], and one-dimensional systems [51]. In such cases, the perturbative expansion often proves inaccurate or inconvenient.
Instead, more symmetric spin representations such as Schwinger bosons [9,52] and slave particle approaches [53,54] are commonly used, but require the introduction of auxiliary fields. Other fermionic approaches include the Jordan-Wigner representation of spin-1/2 operators [55], and its generally complicated-to-use generalizations to higher dimensions [56,57] and higher spin [58][59][60]. An important, equivalent alternative to the HP representation that also uses a single boson species, but avoids the square-root, is the Dyson-Maleev representation [61][62][63][64]. However, it has the drawback of generically breaking hermiticity. This is by no means an exhaustive list of spin representations -indeed, other representations can be be found in Refs. [65][66][67][68][69][70]. Since each of the available approaches has its own unique advantages and drawbacks, we will in this paper derive expressions for spin operators that i) involve only one boson species satisfying the canonical bosonic commutation relation, ii) preserve hermiticity, and iii) do not include square-roots or other non-polynomial functions of operators.
Some of the expressions we derive were previously found to finite order in Ref. [71,72] by a matching matrix elements (MME) method and have also been usefully applied in [73] to capture effects beyond the reach of a 1/S expansion. Unlike the normal Taylor expansion of the HP OSR the MME expansion and our result are able to correctly describe the symmetry in a Heisenberg model with easy-plane anisotropy as we will see later in the text. This is a long standing problem and was discussed using a slightly different approach in [74]. Our expansion thus naturally captures the same physics. However, unlike previous works, we present results to all orders and show that the expressions are exact when truncated to an appropriate order that depends on the spin length S. This feature was missed in all previous discussions we are aware of, since they focused entirely on reproducing commutation relations of spin operators. We, however, use a slightly softer exactness criterion. Namely, we require only that the operators are block-diagonal with physical and unphysical subspace blocks. In the commutator language we require that the commutators are reproduced up to a term that acts exclusively on the unphysical subspace, without coupling to the physical subspace. This is akin to allowing an inaccessible "dark sector" in the spin operator algebra. A more detailed discussion of this rationale is given in the main text.
The manuscript is structured as follows. In the next section of the paper we discuss how to compute squareroots of operators by using a differential equation approach. In section IV we show how this formalism may be used to find a series expansion for √ a † a near a † a ≈ 0 in terms of integer powers of (a † a), which is an unexpected result because √ x cannot be expanded in integer powers of x near x = 0. Of course, since a † a is an operator a † a ≈ 0 is a shorthand for "in the part of the Hilbert space where where matrix elements are close to zero". We will use similar shorthands throughout the text. This shows that a Taylor series may not always be ideal for finding power series expansions of operator functions. In section V we then apply the method to the Klein-Gordon particle in a magnetic field with small or zero-masssuch as in graphene. In section VI we present our main application to the Holstein-Primakoff representation of spin operators. We stress that the results we obtain are exact expressions for spin operators that are polynomial in bosonic operators. Lastly we present our conclusion.

III. GENERAL FORMALISM
The goal of this section is to find an operator differential equation that can be used to calculate a square root of two operators, where O 1 and O 2 are both operators defined on the same complex Hilbert space H. We will make two simplifying assumptions. First, we assume that a square-root of one of the operators, O 1 , is known or easy to calculate. Second, we assume that the two operators commute, [O 1 , O 2 ] = 0. Both these assumptions also have to made in order for a Taylor expansion in O 2 to be viable (the more generic case is more involved, see Appendix A for details). For instance one could have O 1 = c1 with c ∈ C, and O 2 any other operator. It should be noted that the OSR of an operator O 1 can have multiple branches, which may seem like an ambiguity. However, the choice of branch will be encoded in the initial conditions for the differential equations we derive, and should be informed by the problem at hand. Different branch choices can lead to different physicse.g. a branch with complex eigenvalues could not be used to describe a Hermitian Hamiltonian. The way we will compute √ O 1 + O 2 is by introducing the second operator O 2 in infinitesimal steps. To keep track of the steps we introduce a dummy parameter s and define Using the assumption [O 1 , O 2 ] = 0 we find that sending if δs is infinitesimal and we therefore did a Taylor expansion of the right hand side. A Taylor expansion of the left side gives us the differential equation that makes it possible to find O √ (s) by introducing O 2 via infinitesimal steps. Note that this also means that √ O 1 for the branch used needs to be invertible, or at least the limit of an invertible operator as we will see in the upcoming section.
The issue with this equation is that calculating the inverse of an operator is difficult. That is, we cannot easily make an ansatz for O √ (s) = n c n (s)Ô n as a sum of operators with s-dependent coefficients and solve this equation because calculating the inverse of the ansatz is difficult.
This issue can be resolved with a little bit of extra work. We define In this case Eq. (3) becomes and we now need to find a differential equation for O −1 √ (s), which can be obtained by Taylor expanding O −1 √ (s + δs) in a similar way to above, One may insert Eq. (5) in Eq. (6), and we find after rearranging that The equation in this form is now useful to find coefficients C i for an ansatz O √ (s) = n c n (s)Ô n because powers of this operator are trivial to compute.

IV. EXPANDING THE SQUARE-ROOT OF THE NUMBER OPERATOR
We may now use equation (7) to find an expansion of √ a † a. In the language of the previous section for this case O 2 = a † a and O 1 = 0 + 1, where 0 + signifies a dummy variable that eventually will take a directed limit to zero. One can make the ansatz and compare coefficients of (a † ) n a n to find a set of differential equations for C n . If we truncate at third order we find and initial conditions The initial conditions were found by comparison to the infinitesimal case that is accurately described by a first order Taylor series. Note that the term with 1 2 √ 0 + represents a directional limit that has to be taken at the end but 0 + can first can be replaced by a dummy variable.
If we solve the equations and set s = 1 and take the One should note that this expression can be put in terms of powers ofn = a † a and is valid near a † a = 0. More precisely, in what sense does this expansion converge to the correct operator? The answer is that by including terms up to a † n a n the n + 1 lowest eigenvalues are exactly reproduced; higher eigenvalues are approxi-mated more accurately as well.
It is important to stress that the square-root √ x is non-analytic near x = 0. Yet we were able to find an expansion in terms of powers of x that is valid near x = 0. Since x = a † a is an operator with a spectrum that has more structure than a scalar, the flow equations for the couplings do not merely reproduce a Taylor series, but rather lead to an improved approximation. For this specific case, we may intuitively see that the improved expansion is possible because the operator only takes discrete values.

V. APPLICATION TO THE KLEIN-GORDON SQUARE-ROOT
The method for finding a non-perturbative expansion of an operator square-root can of course also be used for the Klein-Gordon square-root Hamiltonian for relativistic particles. Let us for instance consider the 2D Hamiltonian If this system is subjected to a constant magnetic field given by A = B 2 (−y, x) one may introduce the magnetic field by minimal substitution p i → Π i = p i − A i and one can introduce creation and annihilation operators,

to find the Hamiltonian
where we introduced a short-hand S = m 2 +|B| 4|B| . The operator now bears a striking resemblance to the squareroot that appears in the Holstein-Primakoff spin representation, which we will discuss later. A straightforward Taylor expansion of the square-root in terms of 1/S already yields corrections to what one would expect from the non-relativistic relativistic limit of large mass This approximation lifts the restrictions to large masses from the non-relativistic limit as long as one considers strong magnetic fields. However, we can do better without the introduction of further complications. That is, we can make the ansatz √ 1 + sa † a = n C n (s)(a † ) n a n , which means that we can employ the first two differential equations from (9) to approximate the square root. For this we have to choose slightly different initial conditions than previously, C 0 (0) = 1, C 1 (0) = 0, C 1 (0) = 1/2 and let s run up to s = 1/(2S). The initial conditions are again found by comparison to a first order Taylor expansion. The result we find is This new approximation is now more reliable for small |B|, m and level number n. This is seen most easily in the case of V (x, y) = 0 where it is easy to check that it reproduces the lowest two energy levels n = 0, 1 exactly (recall thatn = a † a = SΠ 2 /(2|B|)).
The advantage of this approximation over an exact solution is that a quadratic V (x, y) can be added and an analytic solution of this approximate problem is still possible because this is still a harmonic oscillator. Note that in this case we would be able to find an approximation that is non-perturbative in 1/S.

VI. RESUMMED HOLSTEIN-PRIMAKOFF EXPANSION
We will now turn to our most interesting application -an expansion for the square-root in the Holstein-Primakoff representation of a spin operator.

A. Review of the method
The Holstein-Primakoff representation [8] of spin-S operators is given as A few notes are due. For finite S only finitely many bosonic excitations correspond to physical states. That is, bosonic excitations correspond to spin projections i.e. S z can only take eigenvalues in {−S, −S + 1, . . . , S}.
Hence, for spin S we have the restriction a † a ≤ 2S, which is also signaled by the fact that the square-root becomes imaginary for higher occupation numbers. This means that the Hilbert space is a Fock space, F (H) = ∞ n=0 SH ⊗n , where S is the symmetrization operator, H ⊗n denotes n tensor products of the single particle Hilbert space H. For spin S the physical part of the Hilbert space is restricted such that it has the basis {|0 , ..., |2S }.

B. Exactness of the Holstein Primakoff approximation
To see that the Holstein-Primakoff representation is an exact description of spin operators it is enough to check that it fulfills the correct spin algebra. For instance [S + , S − ] = 2S z . This reasoning is slightly restrictive so let us soften it a bit. The key feature of the Holstein Primakoff representation is that the spin operators S +,−,z reproduce the exact spin operators on the physical part of the Hilbert space and at the same time have no elements that couple to the unphysical part of the Hilbert space. That is, in the occupation basis they have the form In particular, for S = 1/2 the explicit form of S + in the occupation basis is One sees that it splits into the physical (highlighted in blue) and unphysical (red) Hilbert spaces like in equation (18).The physical block is just the conventional S + matrix for spin 1/2. Importantly there is no coupling between physical and unphysical parts of the Hilbert space. This is what makes the method exact.
Note that, because of this block structure, a spin Hamiltonian exactly written in the bosonic language will also separate into physical and unphysical blocks because the product of block diagonal matrices stays block diagonal. That is, the Hamiltonian is block diagonal of the form One now can see that diagonalising the Hamiltonian one will find the exact physical eigenvalues and spurious unphysical ones.

C. Usual Approach: Taylor expansion
While the expressions in Eq. (17) provide an exact way to represent the spin operators this is not too useful by itself because the square-roots are impractical to work with. One usually does a Taylor expansion around large S, using 1/S as expansion parameter.
This approach is most often also used in the case of S = 1 2 , where it is slightly surprising that it is justified. To see why it is justified recall that as mentioned above for smaller spins only states with few bosonic excitations e.g.
{|0 , |1 } for spin 1 2 are physical. Therefore acting in this part of the Hilbert space a † a| phys ≤ 1 and the expansion is valid.
Although the expansion is useful it is not exact when truncated at any finite order. The spin operators S +,− no longer separate into physical and unphysical blocks, but couple physical and unphysical parts of the Hilbert space. For example, for spin 1/2 the spin operator S + in Eq. (21) has the form One may see that physical and unphysical parts of the Hilbert space get coupled by the term 5 8 √

.
Generically a spin Hamiltonian using this approximate bosonic language when expressed in occupation number space has the form where ∆ is the small coupling between physical and unphyical parts of the Hilbert space. It leads to unphysical contributions in the physical eigenvalues. The method is not exact anymore.

D. Improved Expansion
As mentioned we can improve on the expansion. One may use the differential equation (7) to find such an improved expansion of the square-root. Like previously one may use the ansatz (8) to introduce − a † a 2S by infinitesimal steps. However, because we need to decrease terms under the square-root rather than increase them, one has to replace d/ds → −d/ds in (9). The second thing that changes compared to before are two of the initial conditions while the other initial conditions in (10) remain unchanged. The solution to these differential equations (9) for s = 1 S with the new initial conditions gives us an improved Holstein-Primakoff expansion up to third order, which we will not present here.
Rather with additional work one may find that it is possible to construct higher order terms by the same scheme. After analysing additional orders one can see a pattern emerge. We find that the full expansion is given as Q n a † n a n a; Q 0 = 1, where we prove later that this amounts to exact expressions for spin operators.
Let us for now truncate at n max = 1 to find and discuss the case of spin S = 1 2 to most easily see what kind of improvement we achieved. One may note that an expansion around large S gives back the results for the Taylor expansion. In that sense our new expansion is a resummation of the Taylor series.
In the occupation basis we find Therefore the spin operator reproduces the physical matrix elements of S + , and the physical block does not couple to the unphysical block like in Eq. (18). One can show also more explicitly that there are no coupling between physical and unphysical parts of the Hilbert space In the same sense as before this method therefore allows us to reproduce the exact eigenvalues of the Hamiltonian. In this sense it is exact.
Of course, this first truncated expression is not exact for higher spins S because couplings to the non-physical states reappear. We can obtain exact expressions also for S > 1/2 by setting n max = 2S. Similarly to the S = 1/2 case these expressions reproduce all physical matrix elements. The proof is given in the Appendix B, but is essentially the same as for spin 1/2. A list of explicit expressions for spin operators up to S = 3 are given in Appendix C.

E. Commutator properties and exactness for the improved expansion
One may ask what happens to commutators. Here the spin 1/2 case again is instructive, (29) While the commutator is not exactly reproduced we can immediately recognize that this is not important because the extra term a † 2 a 2 does not couple unphysical and physical parts of the Hilbert space and solely affects the unphysical parts. It is therefore of no physical consequence.
This additional term often was understood as rendering the expressions for spin operators approximate [71]. After all, the most commonly used criterion for ruling out if an operator can be expressed in a certain way is by checking the commutation relations. Here we stress that this criterion can be softened. Namely it can be enough to reproduce commutation relations up to the addition of a term that acts solely in the unphysical part of the Hilbert space and does not couple to the physical part of the Hilbert space.
In some cases more stringent exactness criteria have been applied, such as requiring that all non-physical matrix elements vanish [70]. This type of criteria can simplify formal quantum statistical treatments, since one does not have to be careful about excluding non-physical states in sums over states. However, in practice these approaches are cumbersome because the associated expansions are infinite and more complicated. Therefore this is only an advantage at the purely formal level.

F. Additional properties of the expansion and comparison to other expansions
One may wonder how this expansion compares to a more conventional Dyson Maleev expansion with S + = a and S − = a † (2S − a † a). Our method has the advantage that S + and S − are treated on the same footing and therefore are related by conventional Hermitian conjugation. This guarantees that the approach will not break hermiticity in the conventional sense, unlike the Dyson-Maleev expansion.
Next one may wonder if an additional perturbative expansion around classical spin configurations may be stacked on top of the expansion as it is done for the more conventional 1/S expansion in non-linear spinwave theory [45]. One may therefore be tempted to identify δ = 1 − 1 2S − 1 in Eq. (26) as an expansion parameter since it corresponds to fluctuation corrections around a classical ground state. That is, one would write S + ≈ √ 2S 1 + δa † a a. This, however, is not possible and becomes clear if one considers that S z = 1/2[S + , S − ]. Then, one can write (30) We find that the physical part of S z has contributions from different orders of δ. This of course means that any expansion in δ will treat S z and S x,y on unequal footing, even at low orders in such an expansion. This, for instance, will result in an unphysical breaking of symmetries in a Heisenberg model or similar even at the lowest order expansion. Therefore, δ cannot be used as an expansion parameter. Additionally, there is no other obvious choice of expansion parameter and ad-hoc expansions in powers of a also lead to unphysical results in non-linear spin-wave theories. Therefore, it seems that the expansion does not allow for an additional perturbative expansion in terms of fluctuations around a classical spin configuration. A mean field theory treatment must include all the terms that are needed to accurately describe spin S for each operator S + and S − .

G. Symmetries and exact properties in the improved expansion
To study symmetries in the new expansion we consider the Hamiltonian for the Heisenberg model with easyplane single-ion anisotropy, Let us first recognize that for S = 1/2 the singleion anisotropy (S x i ) 2 should result in a trivial number (S x i ) 2 = 1/4 that does not affect the spin-wave excitation spectrum. However, in the usual Taylor expansion with S + ≈ a − 1 2 a † a 2 one finds that which has unphysical contributions in the physical part of the Hilbert space, e.g. a † a. In other words, the Taylor expansion introduces unphysical artifacts. In the new expansion for S = 1/2, however, we have S + = a − a † a 2 and find that (33) The additional non-constant terms we find have non-zero contributions only in the non-physical part of the Hilbert space, and do not couple to the physical part of the Hilbert space. This can easily be verified explicitly by computing the operator in the occupation number basis using Eq. (27). The non-physical terms are therefore of no consequence for physical states, and could just as well be dropped. This means that the new expansion properly reproduces the fact that (S x i ) 2 contributes only a trivial scalar for spin 1/2.
Next we recall that the Hamiltonian is symmetric with respect to the symmetry generated by the generator g = i S x i , i.e. C = [H, g] = 0. Again we will only check spin 1/2 for simplicity, but similar results will hold for higher spins. Let us first see what happens if we use the usual Taylor expansion approach to compute the commutator. We find that where (i) ↔ (i + 1) is a shorthand for the same terms with i and i+1 switched. Here we can see that the operators in the first line couple the physical two site states in } to non-physical two-site states in A np = |n i |m i+1 |(n > 1) ∨ (m > 1) , which are the states where at least one of the two sites is more than single occupied. The symbol ∨ denotes the inclusive "or" (disjunction) operator.
For the new expansion, on the other hand, we find that From here one may observe the term a † j 2 a 2 j to see that only matrix elements for the unphysical states A np will be non-zero. The operator also does not couple to the physical states in A p . We can hence conclude that, unlike the Taylor expansion, the bosonic expressions for the spin operators in the new expansion do not break symmetries present in the original spin operator language.

VII. CONCLUSION
We were able to demonstrate the surprising result that the square-root of an operator Ô may be expanded in an integer power series aroundÔ = 0. We believe that the approach can be usefully applied to other operator square-roots in theoretical physics and that the observation is useful for finding better expansions of other operator functions where a Taylor expansion fails.
The methods described in this paper allowed us to find a significant non-perturbative improvement on the Taylor expansion for the Holstein-Primakoff realization of spin operators. We expect these results to be useful to better treat spin models in different mean field approaches if there is no clear classical spin configuration around which one could expand. We therefore hope that the approach will prove useful for the study of spin liquid phases. ACKNOWLEDGMENTS We thank C. D. Batista and G. Marmorini for useful discussions. M.V. and G.A.F. gratefully acknowledge partial support from the National Science Founda-