Exploiting polarization dependence in two dimensional coherent spectroscopy: examples of Ce 2 Zr 2 O 7 and Nd 2 Zr 2 O 7

,


I. INTRODUCTION
A defining feature of quantum spin liquids is the presence of fractionalized excitations, which cannot be created on their own.When probed with linear response methods, such as inelastic neutron scattering, this leads to signal which is broad in wave-vector and frequency, corresponding to the many ways of distributing the scattered momentum and energy amongst the multiple fractional excitations [1][2][3][4].The lack of a sharp signal leads to a challenge in the experimental study of spin liquids: how can continua of fractionalized excitations be distinguished from broad scattering due to the simple absence of longlived excitations?Two dimensional coherent spectroscopy (2DCS) offers a resolution to this conundrum.2DCS probes the nonlinear response of the system to a pair of radiation pulses, at different times.By varying the temporal separation of the pulses and the measurement, one can measure the response as a function of two time coordinates, or, when Fourier transformed, two independent frequencies: ω t and ω τ .Fractionalized excitations show up as streaks of intensity in the ω t − ω τ plane.Of particular importance is the rephasing signal: a streak along the ω t = −ω τ direction [5].This signal is important because it is robust to broadening effects, and allows such broadening effects to be distinguished from extended signals from fractionalized continua: the broadening due to finite lifetimes appears transverse to the streak, while the intrinsic continuum width is given by the length of the streak.It thus becomes possible to directly separate out fractionalized quasiparticle lifetimes, and identify them as long-lived excitations.Since the original proposal to use 2DCS as a probe of fractional excitations [5], several theoretical calculations have demonstrated its usefulness in principle for various models in 1 and 2 dimensions [6][7][8][9][10].Some of the most promising spin liquid candidates are found in the rare-earth pyrochlore oxides R 2 M 2 O 7 [11].In recent years, Ce-based pyrochlores Ce 2 M 2 O 7 (M=Zr, Sn, Hf) have emerged as particularly likely candidates [12][13][14][15], and much effort has been devoted to understanding their correlations and microscopic coupling parameters [16][17][18][19][20].As with other spin liquid candidates, the question of how to interpret inelastic scattering data, in the absence of sharply defined excitations, is a key issue.
These materials also exhibit interesting responses to externally applied magnetic fields, due to the highly anisotropic nature of the crystal electric field (CEF) environments at the rare-earth sites.In particular a magnetic field aligned with the [110] crystal axis splits the system into two sets of spin chains, α and β, with the α chains being polarized by the applied field while the β chains feel neither the effect of the field, nor the neighbouring α chains.This emergence of effectively one-dimensional physics out of a three-dimensional correlated system has been observed in the quantum pyrochlores Ce 2 Zr 2 O 7 [21] and Nd 2 Zr 2 O 7 [22], as well as in the canonical classical spin ices Dy 2 Ti 2 O 7 [23,24] and Ho 2 Ti 2 O 7 [24,25].
In this work we investigate how variations in the polarization of a probe field can be used to extract different information from the 2DCS response of a system, focusing on the quantum pyrochlores Ce 2 Zr 2 O 7 and Nd 2 Zr 2 O 7 with an applied [110] magnetic field.Their highly anisotropic microscopic Hamiltonians provide a way to explore how the dependence of optical matrix elements on probe field polarization allows one to tune the 2DCS response.The emergence of effectively onedimensional spin chains in an applied field allows us to study this response in a controlled way while remaining in a setting where fractional excitations appear, in this case on the unpolarized β chains.
Further, we focus on providing predictions for 2DCS experiments in the materials Ce  as estimates of the interaction parameters are available [16,19,[26][27][28], and the emergence of chains has been established [21,22].We show how polarization variations allow separate access to the dynamical correlations of both α and β chains, and calculate the dependence of the signal on the nature of the in-field ground state and on the microscopic interaction parameters.For Ce 2 Zr 2 O 7 we find that a [001] polarization of the probe field is particularly powerful for probing the proximity to a nearby quantum critical point, predicted from the parameter values in [16,21].This is because the [001] polarization suppresses the matrix elements to excite all states bar those at the lower band edge of the β-chain spinon continuum.On the other hand, a [1 10] polarization probes the spinon density of states across the full spectrum.We demonstrate that a streak-like response is still observed on the α chains are probed using a [110] field, despite the fact that excitations are not fractionalized on these chains, due to the probe field having no matrix element to excite single magnons.
For Nd 2 Zr 2 O 7 , we find that the 2DCS response is profoundly affected by a strong mixing between dipolar and octupolar degrees of freedom found in the microscopic Hamiltonian [26], expressed through a mixing angle θ.In particular, we find that θ controls the relative contribution of one-and two-magnon states to the 2DCS response of the α chains in a [110] polarization set-up.For the [1 10] polarization, where the experiment probes spinons on the β chains, we find that the relative amplitude of the rephasing signal along ω t = −ω τ is suppressed at values of θ far from a multiple of π/2.
The Article will proceed as follows: in Section II we review some necessary background theory regarding both the materials and the experimental method, as well as presenting some new results regarding the effect of quantum fluctuations on the ground state phase diagram; in Section III we give details of our analysis, which makes use of several methods in various limits of the problem (Jordan-Wigner fermionisation, Hartree-Fock mean field theory and perturbation theory); in Section IV we show the results of our calculations as they apply to Ce 2 Zr 2 O 7 and Nd 2 Zr 2 O 7 , before concluding in Section V.Both Ce 2 Zr 2 O 7 and Nd 2 Zr 2 O 7 are rare-earth pyrochlore magnets of a dipolar-octupolar character.In either material, the magnetic ions (Ce 3+ or Nd 3+ ) lie on a pyrochlore lattice, consisting of corner sharing tetrahedra.The geometry of the pyrochlore lattice is displayed in Fig. 1.There are two kinds of tetrahedra in the lattice, related by inversion symmetry, conventionally referred to as 'A' and 'B' or 'up' and 'down' tetrahedra.Each site is shared between one 'up' and one 'down' tetrahedron.The lattice can be split up into four fcc sublattices L i (i = 0, 1, 2, 3), each generated by taking one of the vertices of a single 'up' tetrahedron and acting on it with the translational symmetries of the lattice.We define vectors ẑi to be unit vectors pointing from the centre of the host 'up' tetrahedron to the given lattice site.
Each rare-earth ion has a fixed total angular momentum J, determined by Hund's rules.The crystal electric field at each site breaks the degeneracy of the 2J + 1 different angular momentum states m J of each ion.Both Ce 3+ and Nd 3+ are Kramers ions, and so the lowest energy states arising from this splitting form a degenerate time-reversal paired doublet.For rare-earth pyrochlore magnets with Kramers ions, this lowest energy doublet can be of two types depending on the transformation properties of the doublet under the double group of point symmetries and time-reversal [11,29].One finds that the doublet states either transform in the same way as spin-1/2 states under the action of the symmetry group, or in the more exotic 'dipolar-octupolar' [29] fashion seen in our materials of study.
In dipolar-octupolar compounds, one can construct a basis of pseudo-spin operators S x , S y , S z acting on the lowest energy doublet such that the magnetisation operator µ projected into this subspace is diagonal, and proportional to S z .
One finds that the S ± in this basis correspond not to physical dipole operators, but to magnetic octupoles.Nevertheless, S x in this pseudospin basis transforms identically to S z under the symmetries of the crystal, and transforms differently from S y [29].Thus the S x pseudospin operator, despite physically being a magnetic octupole operator, is often referred to as 'dipolar'.The degeneracy in the set of crystal field doublets is lifted by interactions between ions on the pyrochlore lattice.Taking into account all symmetry allowed interactions between nearest neighbour pseudospins, one arrives at the following form for the Hamiltonian: Note that we are allowed a coupling between S x and S z due to their identical transformation properties under the symmetries of the crystal.We have also allowed for a uniform magnetic field H, which couples only to the magnetic dipole operator S z .The S z -S x coupling terms are often difficult to work with directly, so we elect to perform the following rotation of our pseudospin basis:: Crucially, in this new basis, both S x and S z have some dipolar character, and so both couple with varying strengths to the external field.S ỹ naturally remains purely octupolar.In the new rotated basis, the Hamiltonian now reads: The angle θ describes the symmetry-allowed mixing of dipolar and octupolar degrees of freedom, and plays an important role in the response of the system to external probes.
In this paper we are interested in phases that exist within Ce 2 Zr 2 O 7 and Nd 2 Zr 2 O 7 when they are placed within a [110] magnetic field.Such a field is orthogonal to ẑ1 and ẑ2 , giving rise to chains of pseudospins within the pyrochlore lattice that couple to the external field (sites on sublattices 0 and 3 which are termed α chains), and chains of spins that do not (β chains).Fig. 1 highlights these chains on the pyrochlore lattice.
We shall here briefly review the types of phases that are predicted by classical mean field theory for the Hamiltonian Eq. ( 6) in the [110] field [30].At zero field, the classical ground state can be in one of six possible phases.Three of these are spin-ice like regimes SI λ (λ = x, y, z), where pseudospins obey a '2-in-2-out' rule with respect to the dominant pseudospin component.For two of these regimes, SI x and SI z , this ordering is of mixed dipoleoctupole character, whilst whilst (SI y ) hosts purely octupolar ordering.The system is in one of these regimes if one of the coupling constants is sufficiently large and positive: The remaining three phases host 'all-in-all-out' type order, and again there are three flavours of this phase depending on the dominant exchange coupling, with AIAO y hosting purely octupolar order.These phases exist when: Fig. 2 and Fig. 3 present phase diagrams for the Hamiltonian of Eq. ( 6) with a [110] field, for model parameters close to those of Ce 2 Zr 2 O 7 and Nd 2 Zr 2 O 7 .It is found within classical mean field theory that the AIAO λ phases remain stable to a finite field, whilst the SI λ regimes are replaced by CHAIN λ phases, wherein the system orders along the α and β chains.Within the CHAIN phases, the classical ground states still obey the '2-in-2-out' rule, but the ground state degeneracy becomes subextensive due to spins on the α chains preferentially orienting themselves parallel to the applied field.With the α chain spins pinned by the external field there are two spin configurations for each β chain that preserve the ice rules on the tetrahedra they intersect.This leaves an Ising degree of freedom on each β chain in the ground state.This is at least the picture for the CHAIN x and CHAIN z phases.The S ỹ pseudospins do not however couple to the external field, and so one finds for low fields that the α chains also retain an Ising degree of freedom in the ground states [30].When the field strength becomes sufficiently large, the coupling between the field and the small x and z components of the pseudospin dominates, and the α chains are again pinned.The first phase at low fields is denoted CHAIN * y , which has a second order phase transition at some critical field to the CHAIN y phase.It is also worth noting that at special values of FIG. 3. Classical mean field ground state phase diagram for parameters close to those of Nd2Zr2O7, as determined in [26], for a finite [110] field.For small to moderate fields, the system retains AIAOz ordering, however a phase transition is found for fields of approximately 0.6 Tesla to the CHAINx phase.We find that, in the regime for which the spin chain model is applicable, the coupling to the field is already comparable in size to the largest coupling constant, Jx.
theta, θ = 0, π/2, π, 3π/2, the coupling between the field and one of S x or S z vanishes, allowing for the possibility of CHAIN * x or CHAIN * z regimes at these points.It is important to note that, at the mean field level, interactions between spins on differing chains are frustrated, and cancel out within each of the CHAIN ground states.One can also demonstrate that the lifting of the ground state degeneracy in the CHAIN phases through quantum fluctuations occurs at sixth order in perturbation theory, resulting in a very small energy splitting [30] much less important than the effects of quantum fluctuations within each chain, which will consider in detail in this paper.for Nd 2 Zr 2 O 7 [26].
Looking first at Ce 2 Zr 2 O 7 , we find that the experimentally estimated parameters lie close to the boundary between the SI x and SI y regimes in zero field.If we enforce strictly that θ = 0, then for finite field the T = 0 classical phase diagram is as presented in Fig. 2(a).The two red lines present the evolution of the two experimental possibilities (permutations of J x and J y ) for the location in parameter space with increasing H.In either case, we evolve through a CHAIN * λ regime until a critical field is reached, at which point there is a second order phase transition into a CHAIN phase with fully polarized α chains.The CHAIN * x regime is however unstable to any finite change in θ, as shown in the second panel of Fig. 2. For θ = 0.01, the system transitions to a CHAIN x phase immediately with increasing H if J x > J y .If J y is instead the dominant positive exchange parameter, then as the dominant interaction is octupolar, we still observe a finite range of h for which the CHAIN * y phase is found, before a transition to the pinned α chain phase.
In the third panel of Fig. 2 we present the phase diagram predicted for θ = 0 if quantum fluctuations on each chain are considered.Further details on these calculations will be presented in Section III B. Using the result that interactions between spins on neighbouring chains cancel out for classical ground states, and that energy splittings between these ground states is introduced at high orders in perturbation theory, we take each chain to be an independent quantum spin chain.The phase of the ground state on the β chains is obviously independent of the applied field, and so the quantum phase transitions observed occur, as a function of h, on the α chains.
We find that the effect of these quantum fluctuations is to drastically lower the critical field at which the transition between the CHAIN * λ and CHAIN λ phases occurs.At J x = J y this line passes through H = 0, indicating that along this phase boundary any finite field will remove the ground state degeneracy of the α chains.
We also see that J x = J y marks a transition between distinct chain phases.This transition is a T = 0 critical point, and is associated to the emergence of gapless excitations on the β chains.The proximity of Ce 2 Zr 2 O 7 to this critical point, and the possible application of 2DCS in probing this is discussed further in Section IV.
Experimental model parameters for Nd 2 Zr 2 O 7 place it within the AIAO z phase, which is stable to finite magnetic fields.For a sufficiently large field however, it has been observed [26], and is predicted in classical mean field theory, that a transition to a CHAIN x phase occurs at finite fields H. Fig. 3 presents the classical mean field phase diagram for Nd 2 Zr 2 O 7 with increasing field strength.Classical mean-field theory predicts a transition at a field strength of around 0.6 T, which is rather higher than found in experiments [22], suggesting that mean field theory overestimates the stability of the AIAO order.We will focus on a high field limit for Nd 2 Zr 2 O 7 , where we are certainly in the CHAIN x phase, and within which we can treat the pinned α chains perturbatively.
For both Ce 2 Zr 2 O 7 and Nd 2 Zr 2 O 7 , placed within a sufficiently strong [110] magnetic field, we then expect to observe chain-like structures, which behave which behave approximately independently from one another on energy scales greater than ∼ 10 −4 J [30], where J is of the scale of the coupling parameters.Within this picture, the leading order quantum fluctuations are excitations which propagate along individual α and β chains, and are either domain walls (spinons), or spin flips (magnons) depending on the quantum phase of that chain.

B. Two dimensional coherent spectroscopy (2DCS)
In systems with fractionalized excitations, or in other situations where particles are created in groups, any zero momentum probe will produce a continuum of different excitations due to the many different ways of satisfying momentum conservation.For the example of pairs of spinons on a one-dimensional spin chain, we find absorption peaks at frequencies ϵ = λ −k + λ k , where λ k is the energy of a spinon with momentum k.In linear response, we thus find a broad continuum, potentially difficult to distinguish from a single peak broadened by life-time or disorder effects.The key application of two-dimensional coherent spectroscopy is to produce sharp signatures of such excitation continua and remove this ambiguity.
2DCS is a non-linear technique which produces response functions of two measured frequencies ω t and ω τ .The desired sharp signatures of continuum excitations are diagonal streaks within this ω t − ω τ plane.This includes a "rephasing" or "spinon echo" streak along the line ω t = −ω τ , for which lifetime effects broaden the response transverse to the peak, while the continuum bandwidth is measured by the length of the streak [5].This allows a clean separation between fractionalisation and lifetime effects.
The procedure we consider is as set out in [5].We couple to the magnetic dipole moments in the magnet through weak, uniform, pulses of magnetic field, and then measure the response of the total magnetisation response to these probe fields along a certain axis.Two pulses are applied, a pulse A at time zero, and a pulse B at time τ .The total magnetisation is then measured at a later time t + τ , denoted M AB (t + τ ).The signature we are interested lies in the non-linear part of the response, thus the linear part of the response must be removed.The magnetisation response to the pulses A and B alone must then be measured and the results subtracted off the full response.If the time integrated magnetic field strengths, in units of ℏ, for the two pulses are A 0 and A τ and we assume that each pulse is sufficiently brief that they can be approximated as delta functions, one finds that the non-linear response of the magnetisation is given by [5]: If the probe fields are assumed to be perturbatively weak, we can evaluate the higher order magnetisationmagnetisation susceptibilities χ (n) using Kubo formulae.
If we measure the magnetisation in the same direction as the probe fields, then we can express these susceptibilities in terms of a single magnetisation M parallel to an as yet unspecified axis.We then define the second and third order susceptibilities as: and Here V is the total number of sites, Θ(t) is a step function ensuring causality and the expectation value is taken with respect to the ground state.Finite temperature effects can be included by replacing this expectation value with a trace of the commutator with a thermal density matrix, which can be evaluated in certain cases using Keldysh contour techniques [9].We find that for the analytically tractable spin-chain models we consider, finite temperatures generically lead to a suppression of the nonlinear response, with the contribution from absorption at a given energy ϵ falling away with T as tanh ϵ 2T .As we are chiefly interested in polarization effects on the shape of the response in the two dimensional frequency plane, we elect to work at T = 0 throughout this work.
The choice of polarization for the probe fields and the measurement axis determines the operator M entering into the Kubo formulae.In the magnetic dipole limit of coupling to an external magnetic field, this will in general be a linear combination of spin operators, and so altering the polarization axis will have highly non-trivial effects on the commutators in the Kubo formulae.In the particular example of dipolar-octupolar pyrochlores, the coupling of magnetic dipoles to a probe field field m is relatively simple; as shown in Eq. ( 1) it is just proportional to m • ẑi S z , where ẑi is the local z axis for each of the four fcc sublattices.The effect of the polarization on the optical matrix elements thus only enters as a numerical factor and the operator content of M is unchanged, which greatly simplifies our analysis.
Even with this simplification, the variation of m •ẑ i has a significant impact on how the 2DCS response probes the excitations of dipolar-octupolar pyrochlores.If we align both of the probe fields, and the axis along which we measure the magnetisation response, with the [1 10] direction then we can couple exclusively to the β chains, whilst a [110] polarization singles out the α chains.
For both of these polarizations, the sign of m • ẑi varies from one site to the next along the relevant chains, however if we choose to perform the experiment with a [001] polarization, we find uniform strength coupling to all sites with constant sign along each α and β chain (but a relative sign between types of chain).Comparing the magnetisation operator in the [001] case to the [110] and [1 10] cases, there is an effective momentum shift of π, and so one can use the [001] polarization to probe different excitation processes within the system.Throughout this work, our analysis of the 2DCS response of our chosen materials will be within an independent chain approximation, and we will focus on the [001], [110], and [1 10] polarization set ups.We will then need the second and third order magnetic susceptibilities of both the α and β chains to a probe field which alternates in sign from site to site: β,alt , and the susceptibilities for a field with uniform sign along each chain: To be precise, for the [110] polarization, the relevant susceptibilities χ (n) ∥ , and the required magnetisation op-erator are defined as: The 'parallel' subscript is used to denote this polarization mode as it is parallel to the constant, uniform applied field H.
For the [1 10] polarization, we similarly define χ ⊥ and their relevant magnetisation operator to be: Finally, for the [001] polarization, the susceptibilities χ (n) z involve contribution from both α and β chains.Both chains have equal magnetisation operators, up to a sign which we factor out of the β susceptibilities: 16) In the following sections, we examine how the 2DCS response varies between these three polarization modes, and how this can be used to probe the properties of our example materials Ce 2 Zr 2 O 7 and Nd 2 Zr 2 O 7 .In the following sections, we shall analyse the functions χ α,direct , and χ (n) β,direct .This analysis will primarily be through a mean field Hartree-Fock treatment, however we make use of a high field perturbation expansion for the α chains of Nd 2 Zr 2 O 7 , as the non-zero value of the dipolar-octupolar mixing angle θ introduces interaction terms to the Hamiltonian that are otherwise difficult to handle analytically.

III. ANALYSIS
A. Non-interacting fermions (Jz = 0, θ = 0) The Hamiltonian of Eq. ( 6) can be mapped to a system of non-interacting fermions for J z = 0, θ = 0 [31], and fortunately this limit is very close to the model parameters predicted for Ce 2 Zr 2 O 7 .Further, the most troublesome non-integrable term for general parameters is the effectively twisted field acting on the pseudospins for non-zero values of dipolar-octupolar mixing angle θ, which does not appear in the Hamiltonian of the β chains.We thus expect the behaviour on β chains at least to also be well described by the exactly solvable limit more generally, provided J z is not too large.
We thus begin our calculation of the required susceptibilities by evaluating them within this exactly soluble limit.This will entail finding the fermionic modes that diagonalize Eq. ( 6) in this limit, and much of this work follows that done elsewhere [5].
We choose always to work in a local spin basis for which the total magnetisation operator M µ,alt or M µ,direct is uniform.In the case of M µ,alt , this requires a rotation of π about the local y axis on every alternate spin site, which both removes the alternating (−1) j factor from the Hamiltonian, and introduces a minus sign in front of J x (J z = 0 for now, but this coupling also has its sign flipped).For the 'alternating' case then, the exactly solvable Hamiltonian expressed in this new basis reads: We have taken out a factor of 1/2 from each pseudospin operator, and have retained the definitions of the parameters in Eq. ( 6) for clarity.The magnetisation operator in this new basis is: To evaluate χ β,alt , we need to time evolve Eq. ( 19) and evaluate the nested commutators in the Kubo formulae.For Eq. ( 18), the fermionization procedure follows that given in [5], which we will only sketch here.One can map Eq. ( 18) to a quadratic fermion Hamiltonian via the Jordan-Wigner transformation.This is then diagonalised using a suitable Bogoliubov transformation, with Bogoliubov fermions defined in terms of the Jordan-Wigner fermions as The Hamiltonian is reduced to: with the functions sin ϕ k , cos ϕ k , and λ k given by The time dependent probe field operator M (t) can be expressed in terms of the c k fermions in the same way as outlined in [5], and then inserted into either Eq. ( 10) or (11) to determine the non-linear susceptibilities.Dephasing/depolarising and lifetime effects can also be included phenomenologically by using a quantum channel model for time evolution, as described in the supplementary material for [5].Including these effects, one arrives at the following expressions for the second and third order susceptibilities in the limit J z = 0, θ = 0: and with Here T 1 and T 2 are phenomenological time scales characterising depolarization and dephasing respectively, and are as defined in [5].At this exactly soluble point, we note that the second order susceptibilities vanish if h = 0, and so χ β,alt always vanishes, since the beta chains don't feel the field.
We note here that the Jordan-Wigner fermion construction is valid in both low and high field phases, and so equally well describes spinons for small h and magnons for large h.However, only spinons are fractionalized excitations, whilst continuum responses are predicted in all limits.The existence of the continuum only indicates the creation of multiple particles in the interaction with the probe field, and does not necessarily rule out singly created particle excitations.As discussed in [5], excitations of single particles can be created in high field limits of quantum spin chains by probing the system using a transverse field (S x or S y ).For dipolar-octupolar pyrochlores, such processes would occur through magnetic octupole coupling, or magnetic dipole coupling for non-zero mixing angle θ, as discussed in Section III C. Also note that in the high field limit the function sin ϕ k falls away as J µ /h, and so all parts of the response diminish at least as fast as J 2 µ /h 2 , with A k and A k falling off with J 4 µ /h 4 .The entire response formally vanishes in the infinite field limit, where excitations are non-dynamic magnons (single spin flips).
We now move to consider χ , for which we effectively have an alternating external field h but a uniform probe field along each chain.χ involves no external field, so can be calculated directly using the above results with the sign of J x reversed.Accordingly, χ β,direct also vanishes.For χ (n) α,direct however, we can either work with an alternating field in the Hamiltonian, or put the alternating sign in the magnetisation operator.We choose to do the former, which will result in a two band formulation of the dispersion.The fermionized Hamiltonian in k-space is found to be: where ϵ k = (J x + J y )/2 cos(k) and ∆ k = (J x − J y )/2 sin(k).One finds that the alternating field term in the Hamiltonian introduces couplings between k and k − π momentum fermion modes.
To diagonalize this Hamiltonian, one must perform a 4 × 4 Bogoliubov transformation.The required unitary transformation is of the form: with The final Hamiltonian is then: The new fermions c k are only defined for momenta in the range k ∈ [−π/2, π/2].The Bogoliubov rotation angle, and fermion dispersion relations are found to be: As this is ultimately the same α chain Hamiltonian diagonalized in a different basis, the two band spectrum is entirely equivalent to the single band case.
Translating the magnetisation operator M α,direct into this fermion basis, one obtains: where we have introduced the notation φk = (ϕ k,+ + ϕ k,− )/2, and The magnetisation operator now creates fermions with equal and opposite momenta from different energy bands, for a total excitation energy of 2 λk = λ k,+ + λ k,− ,and can also facilitate transitions between bands, with and associated energy change ∆λ k = λ k,+ − λ k,− .Thus, when the second and third order susceptibilities are recalculated, we pick up new time dependencies arising from inter-band transitions.We also find that the second order response vanishes, just as it does in the h = 0 limit, leaving one only with the third order contribution: A (1) = − 16 sin 2 φk cos 2 φk sin( 2λ Here we introduce the notation φk = (ϕ k,+ + ϕ k,− )/2.When h = 0, the two bands are exactly degenerate (equivalent to the symmetry of the positive momentum spectrum about k = π/2 in the rotated basis), and so we recover the zero field result.Taking T 1 and T 2 to be phenomenological constants independent of other model parameters, we then insert the same exponential in time factors as in Eq. ( 24) and Eq. ( 25), so as to be consistent with the zero field limit.
Whilst the form of the response to a direct rather than an alternating probe field is somewhat altered, it retains the crucial features that make the 2DCS response useful in identifying fractionalized excitations, in particular the rephasing part of the response, A (4) is broadly unchanged, save for the values of the excitation energies.We also note that the additional minus sign in front of J x in χ µ,direct at an exactly solvable limit of the model, with θ = 0 and J z = 0.The latter of these can be relaxed for sufficiently small values of J z , at the expense of a density-density interaction between Jordan-Wigner fermions.If indeed J z is not too large, this interaction can be treated in the Hartree-Fock approximation.This is the method that was used to calculate the phase diagram in Fig. 2(c), and is also used to calculate the 2DCS response of Ce 2 Zr 2 O 7 presented in Section IV.
We choose to work with a basis such that the Hamiltonian takes the form of Eq. (18).Including finite J z , the modified Hamiltonian now reads: We then take the following mean field approximation for the density-density interaction term: (44) Note that one cannot set the anomalous expectation values to zero, as the Bogoliubov rotation angle is non-zero.
The above expectation values are then evaluated using the ground state of the bare J z = 0 model, which has the effect of re-scaling the parameters J x , J y , and h.A new ground state can then be determined, and a new estimate of the expectation values calculated.This procedure then defines a set of iterative equations for our model parameters, and the functions cos ϕ k and sin ϕ k which depend on them: If these sequences converge, then the new model parameters are their limiting values.We refer to the renormalized parameters generated from the Hartree-Fock procedure as J ′ α , h ′ .We importantly find that k cos ϕ k vanishes when h (0) = 0 (equivalent to a half filling of Jordan-Wigner fermion modes) and as a result h = 0 is a fixed point of the iterative equations.No effective field is generated by interactions, however they can re-scale an already existing field.
Thus at this order of approximation, the effect of finite J z is only to re-scale the model parameters.No new processes are introduced, and the 2DCS response is not fundamentally altered.
Further corrections beyond Hartree-Fock can be considered perturbatively, and as shown in [9] these lead to finite fermion lifetimes, and other anomalous broadening effects.For our purposes, finite lifetime effects and other sources of broadening will continue to be considered phenomenologically.
Fig. 4 shows how the model parameters flow for increasing values of J z , both in the absence of an applied field (relevant for β chains) and for h = 1.The red dotted lines show points where a quantum phase transition occurs in the chain model, accompanied with a gap closing in the excitation spectrum.The phase diagrams shown pertain to the individual spin chains, and a much wider range of parameters is presented than those that correspond to the CHAIN phases of dipolar-octupolar pyrochlores.
If h = 0, the quantum spin chains host four distinct phases, corresponding to antiferromagnetic, or ferromagnetic ordering parallel to either the x or y pseudo-spins.For example, for large positive J y , the chains are in an antiferromagnetic-y phase.Note that this notion of ferromagnetic vs. antiferromagnetic ordering is basis dependent.For finite fields, a fifth field polarized phase exists across the J x = J y diagonal, with ordering parallel to the applied field.
We observe that finite J z does not lead to any phase changes if h = 0, provided J z remains small enough for Hartree-Fock to converge.If h ̸ = 0, then the field polarized phase is stabilised for positive J z , and destabilised for negative J z .Further, unlike in the h = 0 limit, finite J z can lead to gap closings, and a transition into or out of the field polarized phase.
For completeness, we also present a formulation of the Hartree-Fock procedure in the basis with an alternating applied field: As the Hamiltonian in this basis is just a rewriting of the constant h field Hamiltonian, the resulting phase diagram and parameter flows are identical to that in the right panel of Fig. 4. Note that whilst we have moved between bases where the Zeeman field is alternating and constant, we have not at any point re-defined our parameters, and so Fig. 4 is the same for both calculations without the need for any further flips in the signs of J x and J z .
C. Effect of θ ̸ = 0: Matrix Pfaffians and perturbation theory The approximation θ = 0 is suitable for Ce 2 Zr 2 O 7 , but Nd 2 Zr 2 O 7 has a significant dipolar-octupolar mixing angle, which cannot be ignored.For the calculation of the β chain responses χ  theta will only be to modify the magnetisation operator M , such that it acts on each site as a mixture of S z j and S x j .We have given an analytic treatment of the susceptibility for purely S z j polarized probe fields, but more general susceptibilities can be constructed in terms of Majorana operators and matrix Pfaffians [5].In this section we shall briefly review this modification to the calculation, before discussing the effects of finite θ on the α chains through perturbation theory in the high field limit.The magnetisation operator acting on the β chains, with a finite mixing angle θ, can be written in a basis such that it takes the form: In general we will find expectation values involving both S z j and S x j operators contributing to the second and third order susceptibilities.However, as the Hamiltonian itself is unchanged, we still retain rotational symmetry by π radians about any of the three principal pseudospin axes, which forces all expectation values appearing in the second order susceptibility to vanish.Using the notation χ (3) zzzz to represent the contribution to the full susceptibility from terms with four S z j operators, and similar for other combinations of z and x operators, we find that χ (3) is now given by: χ (3) = cos 4 θχ (3)  zzzz + sin 4 θχ (3)  xxxx + sin 2 θ cos 2 θ χ (3)  zzxx + χ (3)  zxzx + χ (3)   zxxz +χ (3)  xzzx + χ (3)  xzxz + χ (3)  xxzz .(52) From symmetry considerations, only contributions with an even number of x and z operators can be non-zero.From symmetry considerations, only contributions with an even number of x and z operators can be non-zero.Following the calculation given in the supplementary material in [5], one can express σ x j and σ y j in terms of the Majorana fermion modes α j and β j , defined as: These Majorana fermions obey the following anticommutation relations: {α i , α j } = 2δ i,j ,{β i , β j } = 2δ i,j , and {α i , β j } = 0.In terms of these operators, σ x j becomes: As detailed in appendix A and in [5], the Hamiltonian can be recast in terms of these majorana modes, and the S x and S y expectation values can be evaluated using a Wick expansion.This results in a series of matrix pfaffians involving bilinear majorana expectation values such as ⟨α i (t)iβ j (0)⟩.Evaluating each pfaffian required to exactly calculate non-linear susceptibilities quickly becomes computationally expensive, and so certain approximations have been made when evaluating susceptibilities using this method in section IV, as further discussed in appendix A.
2. θ ̸ = 0 calculation for α chains: Perturbation Theory Moving now to the situation on the α chains, things are somewhat more complicated.The term −h sin θ j S x j is highly non-local in Jordan-Wigner fermions, precluding an interpretation in terms of a simple interaction between particles.Instead, we focus on the high field limit in the original pseudospin basis defined in terms of J xx , J yy , J zz , and J xz [Eq.( 2)], wherein the unperturbed part of the Hamiltonian takes the simple form: The ground state is then simply the state with all pseudospins in the σ z = 1 eigenstate.Periodic boundary conditions are imposed.Excitations from the fully polarized ground state are spin-waves, or magnons, with the magnon excitation energy being h.To account for small but non-zero values for the other coupling constants, we turn to leading order degenerate perturbation theory.
First examining the single magnon sector, the effective Hamiltonian has components where the roman indices i, j denote the location of the single flipped spin (magnon).The eigenstates of the Hamiltonian are simply momentum eigenstates |k⟩.To first order, the energy of the ground state is E Ø = − hL 2 + JzzL 4 , and we find that the magnon excitation energies above the ground state are given by: We can similarly write down an effective Hamiltonian for the two magnon sector.Even in this perturbative limit the magnons are interacting, both through a hard-core constraint, and a non-zero J zz , which gives an energy cost or benefit to domain walls depending on its sign.We find that J zz promotes the existence of eigenstates where the two magnon are primarily on neighbouring sites, and propagate as a single free particle, i.e. as bound magnon pairs.
Speaking qualitatively, the first order correction to the ground state |δØ⟩ includes both unbound zero momentum magnon pairs, |k, −k⟩, and the zero momentum bound state |0⟩ 2,bound .Due to the effect of interactions, the exact first order correction is a little more complex, but the above description is approximately correct.In the limit where J zz = 0, we find that the ground state couples only to the unbound magnons, and that this coupling is proportional to 2 sin(k).The sine function here reflects the fact that the hard-core magnons must avoid being on top of one-another, and so behave in this way like fermions.
Writing a general two-magnon sector eigenstate as |n⟩ 2 , whose energy is ϵ n;2 , and introducing coupling parameters f n = j ⟨j, j + 1|n⟩ 2 (|j, j + 1⟩ is the state with both spins j and j + 1 flipped relative to the ground state), we find that the first order correction to the ground state is: We now proceed in calculating the required non-linear susceptibilities via a Lehmann expansion of the required time-dependent expectation values that arise out of the nested commutators in the Kubo formulae.This entails inserting several resolutions of the identity in terms of energy into Eqs.(10) and (11), such that all time evolution operators act on energy eigenstates.We can then take a time-independent matrix element out of the time-dependent commutator.Let the remaining time and energy dependent factor be written as a function G({t i }, {E i }), then the third order susceptibility, for example, becomes: where Σ z = j σ z j .The calculation then proceeds by expanding each of the energy eigenstates appearing in Eq. (59) in powers of J µ /h.At zeroth order, the Σ z operators act trivially on every eigenstate, and each energy E p , E q , E r introduced in the Lehmann expansion is set equal to the ground state energy.As a result, G is found to be timeindependent and real, and the whole response vanishes.We find at first order in J µ /h that the matrix element is proportional to ⟨Ø|δØ⟩, which is zero, and so again the response vanishes.The first non-zero contribution thus comes at second order.We find that two types of terms are generated, those which involve excitations of k = 0 magnons, and those that excite two magnon states with zero total momentum.Both of these have identical forms for their time dependencies, and we find no rephasing behaviour to leading order in perturbation theory.The full third order susceptibility for constant probe and external fields is found to be: − sin(ϵ 2;n (t 1 + t 2 + t 3 )) .(60) In the limit J zz = 0, f n = 2 sin(k) for all non vanishing f n , and the corresponding ϵ 2,n are approximately equal to 2λ k .This matches precisely the large h limit of Eqs.[( 21)-( 23)] to second order in J µ /h.
We can perform a similar analysis for the case of second order susceptibilities.Again we find that the first nonvanishing terms arise at second order in J/h.In this case, we obtain: The second order response then also has a part due to the excitation of k = 0 magnons, and a part due to the two magnon sector.
Thus the key result for non-zero θ is the possibility of exciting single magnons, rather than pairs of excited particles, demonstrating that the polarized regime does not possess fractionalized excitations.As mentioned previously, single magnon excitations should be observable even with θ = 0 if one couples directly to the pseudospins S x or S ỹ directly.Such processes should be seen for our magnetic probe field if we extend our treatment to include octupolar couplings, however the resulting responses will be much weaker than the dipolar contribution.
We are also interested in the response of Nd 2 Zr 2 O 7 to a [001] probe field.The full response will have contributions from χ α,direct , and χ (3) β,direct .The β chain response can be treated within the Hartree-Fock approximation using the Pfaffian expansion as detailed above and in appendix A. For the α chain responses in the limit of large h, we can make use of a similar perturbative expansion to as above.
The only change to the analysis is to replace the operator Σ z with Σz = j (−1) j σ z j , as we choose to remain in a basis in which z i • H is constant, which causes the direct probe field to acquire this alternating sign.Σz is no longer proportional to the unperturbed Hamiltonian H 0 , and so its action on the unperturbed states is in principle non-trivial.However, if we consider an system with an even number of sites, Σz annihilates the ground state |Ø⟩, and the states with two neighbouring spin flips |j, j + 1⟩.As a result, the only surviving contribution to the non-linear response of the α chains in this limit, to second order in (J µ /h), is the single magnon contribution to χ (3) α,direct .Acting on the |k = 0⟩ magnon state with Σz returns −2 |k = π⟩, and similarly Σz |k = π⟩ = −2 |k = 0⟩.There is then only one non-vanishing term in the Lehman expansion of χ (3) , and we find that the full α chain response to second order in perturbation theory is given by: Here λ 0 and λ π are the excitation energies of k = 0 and k = π magnons respectively.The time dependence of Eq. ( 62) does not depend on only one energy.This is similar to the predicted [001] response of α chains for zero mixing angle given in Eq. ( 38), and again we interpret the modified time dependence as stemming from transitions between magnon modes induced by the effective π momentum probe field.
One thing to note is that we have no contribution at second order to the response from the two magnon sector.Comparing this result to Eq. ( 38), we find that if the high field limit is taken for the θ = 0 exact expressions, the leading order contribution occurs at fourth order in J µ /h, which is consistent with the absence of such terms in this perturbative calculation.

IV. RESULTS FOR Ce2Zr2O7 AND Nd2Zr2O7
In this section, we present predictions for the 2DCS responses detailed in Section III for the two rare earth pyrochlore magnets Ce 2 Zr 2 O 7 and Nd 2 Zr 2 O 7 in the presence of a sufficiently strong [110] field such that mean field theory predicts a CHAIN ground state.We examine the responses of both materials for the [110], [1 10], and [001] polarization channels, and discuss the qualitatively different information provided by these different experiments about excitations in both systems.
Where the Hamiltonian of each α or β chain is tractable in the Hartree-Fock approximation, and θ = 0, analytical results are presented.Otherwise, a mixture of numerical results calculated using the matrix Pfaffian or perturbative methods discussed in Section III C are employed.Note that while the matrix Pfaffian method is in principle numerically exact, we employ several approximations to the exact calculation to improve computation time, as described in Appendix A. These approximations give rise to a discrepancy in the global scale of the response of a factor of order one when compared to the exact calculation, but the functional form of the response in ω t and ω τ is unchanged.Thus the arbitrary units scales of response functions presented in this section may differ between plots where different computation methods have been used.The scale of different panels within the same figure are however always directly comparable.
That we have θ ≈ 0 means that the most troubling nonintegrable term in the Hamiltonian of the α chains, the coupling to S x, is removed and the response from both chains, for direct and alternating probe fields, can be treated within the Hartree-Fock approximation.The rescaled β parameters are independent of the applied field, and are given by: The re-scaling of the α chain parameters is field dependent, as shown in Fig. 6.
For both α and β chains, both the rescaled J x and J y are moderately larger than their bare values, whilst their difference remains approximately constant.As H is increased, the values of J x and J y smoothly decrease, but remain higher than their bare values.The re-scaling of the effective magnetic field acting on each α chain is the most significant modification, with an increase in the effective field strength of approximately 0.06 Tesla observed for an initial field strength of 0.2 Tesla.The Hartree-Fock rescalings are independent of which of J x or J y is assigned the larger value, and so too is the 2DCS response.For definiteness, we take J x = 0.064meV, and J y = 0.063meV throughout.
Before discussing the results in detail, it is worth describing the types of signatures observed in the Fourier transforms of the magnetic susceptibilities χ (3) (t, t, t + τ ) and χ (3) (t, t + τ, t + τ ).Their Fourier transforms are functions of two variables, ω t and ω τ .The signature we are interested in lies in the absorptive part of the response, so we take the imaginary part of χ (3) (ω t , ω τ ).Looking at one term in particular, consider the A (1) term [Eq.( 26)] of χ (3) (t, t, t + τ ) in the limit 1/T 1 = 1/T 2 = 0.This contribution has the following form as a function of time: Fourier transforming this and taking the imaginary part, one obtains peaks in the response at ω t = ω τ = ±2λ k .However, these peaks are not pure delta functions, as one is taking the imaginary part of a product of two complex functions.Generically, the full Fourier transform at a given peak has the functional form: Making use of the identity 1 x−a = −iπδ(x − a) + P 1 x−a , where P denotes the principal value of a function, we find that the imaginary part of a given peak function has the form: The delta function part is formally diverging, and negative, whilst the principal value part of the response peak changes in sign as (ω t − ϵ) or (ω τ − ϵ) swap sign.The A (1) term gives rise to many such peaks along the line ω t = ω τ in the ω t − ω τ plane.For small chains and only small amounts of broadening, each set of peak remains distinct, however with greater broadening and larger system size, the regions of different sign begin to overlap.
Referring to Eq. (65) and panel (a) in Fig. 7, if we follow the ω t = ω τ line, we see that as we cross through each peak generated by A (1) , that the function changes from large and positive, to large and negative at the Lorentzian (broadened delta function) part of the peak, to large and positive again.If the broadening is further increased, we find that this sign variation tends to wash out the signal leaving a featureless broad response, as shown in Fig. 7 (c).
On the other hand, the A (4) term of χ (3) (t, t, t + τ ) [Eq. 29] produces identical peaks, but along the ω t = −ω τ line.Moving along this line, the second term of Eq. (65) is negative, the same sign as the delta functions, and so as one moves along ω t = −ω τ , the sign of the function does not vary, as shown in Fig. 7 (b).The signal resembles a negative 'trench' flanked by positive ridges.It is found that when the broadening on the response peaks is increased, this orientation of the streak is resistant to being washed out in the same way as the A (1) signal, as shown in Fig. 7 (d).This result underlies the focus on this rephasing part of 2DCS responses.

[1 10] polarization: probing β-chain spinon density of states
A particularly interesting property of the Ce 2 Zr 2 O 7 model parameters is how close they lie to the predicted phase transition between the CHAIN x and CHAIN y phases.In terms of the behaviour on the α and β chains, for all but the very smallest fields, the α chains are in a field polarized phase.The β chains are however close to a critical phase transition, with a predicted gap to excitations of approximately 1µeV, corresponding to the excitation of two k = 0 spinons.
Being so close to criticality means that a [1 10] polarized probe field (m ⊥ ), will provide a particularly clear picture of the full spectrum of spinon pair excitations, as sin ϕ k is very close to one for all k.The response to m ⊥ is determined entirely by χ β,alt because the second order response vanishes on the β chains.Referring to Eq. ( 9), two different time orderings contribute to the full nonlinear response, χ β,alt (t, t, t + τ ) and χ (3) β,alt (t, t + τ, t + τ ), both of which are given by the expression in Eq. ( 25), with the functions cos ϕ k , sin ϕ k , and λ k as given in Eq. (21)(22)(23) with h set to zero.
One finds that precisely at the limit J x = J y that cos ϕ k =0 and sin 2 ϕ k = 1 (see Eqs. [( 21)-( 22)]).This FIG. 8.A cut along ωt = −ωτ in the 2DCS response of Ce2Zr2O7.The solid blue curve is the calculated non-linear response, whilst the red dot-dashed line is the density of states for the reported model parameters.A large dephasing time of T2 = 750meV −1 is used to minimise the effects of broadening on the shape of the response, whilst maintaining a smooth curve.We set A0 : Aτ = 1 : 100 (see Eq. ( 9)) so as to single out the rephasing signature.We observe close agreement between the two curves, but note some deviations at the upper and lower band edges arising from the still finite broadening.
means that the rephasing signature would have uniform strength across the full spectrum.In the continuum limit, and for minimal broadening, the variation in the height of the response is approximately equal to the density of states in energy, as shown in Fig. 8. Thus the 2DCS response with the [1 10] polarization provides a close approximation of the spinon density of states on the β chains.
In Fig. 9, we present plots of the response of the β chains of Ce 2 Zr 2 O 7 to a [1 10] probe field.Fig. 9(a) shows the response for the estimated model parameters [16], whilst (b) and (c) demonstrate how the response would vary with an increasing anisotropy J x − J y .The strong diagonal signatures in the lower right and upper left quadrants are the aforementioned rephasing signatures, which run from close to the origin out to ω t = −ω τ = | Jx + Jy |.The increase in intensity that is observed at the upper end of the spectrum is due to the diverging density of states at the upper band edge, as the energy of a spinon pair 2λ k is stationary at this maximum as a function of k.
The horizontal bar observed along the ω τ = 0 axis originates in the χ (3) (t, t + τ, t + τ ) part of the non-linear response.The two different time orderings contribute with different weights, with χ (3) (t, t + τ, t + τ ) contributing with a strength A 2 0 A τ , and χ (3) (t, t, t + τ ) with A 0 A 2 τ , as seen in Eq. ( 9).Thus, in principle, the contribution from χ (3) (t, t, t + τ ), which contains the rephasing term, can be favoured by tuning the pulse strengths [9], provided that the second order contribution also vanishes (as it does for the β chains).
With increasing anisotropy J x − J y , we can see in the 2DCS response that the gap to excitations widens and the bandwidth of pair excitations narrows as we tune away from the critical point.In Fig. 10, we again tune away from the predicted model parameters for Ce 2 Zr 2 O 7 , but this time by introducing a finite mixing angle θ.To calculate the non-linear response, we must make use of the matrix Pfaffian method discussed in Section III C, which is somewhat more expensive computationally, and so we are limited to chain lengths of L = 40, and a coarser frequency resolution.Following a similar procedure as that detailed in the supplementary material of [5], we impose open boundary conditions, and impose a cut off of 5 sites from the end of each chain in sums over sites to avoid the effects of edge states.
With finite mixing angle θ introduced to the Ce 2 Zr 2 O 7 model parameters, we observe that the response becomes dominated by a signal along ω t = 0, whilst the rephasing contribution appears to remain approximately constant in amplitude with changing θ.This is shown in Fig. 10, for mixing angles of θ = 0, θ = 0.3, and θ = 0.5.As the response functions in Fig. 10 come only from the β chains, whose Hamiltonian has no dependence on θ, the effects for finite θ must originate in the changing matrix elements between the probe field and the pseudospin degrees of freedom.

[110] polarization: two-magnon excitations on α chains
If we instead choose to align our probe field and measure the magnetisation parallel to the [110] direction, we can extract the response of the α chains alone.The quantum phase transition on the α chains between a magnetically ordered phase and a field polarized phase is predicted occur at | Jx − Jy | = 2 h, which corresponds to a field of the order H c ≈ 0.001T for Ce 2 Zr 2 O 7 .Thus for any moderately large field, the α chains are field polarized, and relatively far from a quantum phase transition.As these chains are in a field polarized phase, and not one hosting one-dimensional ferromagnetic/antiferromagnetic ordering like on the β chains, the fundamental excitations are magnons, which are not fractionalized.However, for θ = 0, the probe field has no matrix element to excite single magnons in the dipole approximation, and thus we still see only responses from excitations of pairs of particles.Thus the 2DCS response will still produce streak-like continua when probing the α chains, despite the absence of fractionalisation.
In Fig. 11, we present the 2DCS response of the α chains of Ce 2 Zr 2 O 7 to a [110] field for increasing field strengths.We observe that as the field strength is increased from 0.1T, to 0.3T and 0.5T that the two magnon bandwidth narrows considerably relative to the gap, with the 2DCS response at 0.5T collapsing almost to single peaks as the energy scale for exciting two magnons (h) becomes much larger than the energy scale governing their dynamics (J µ ).25), within the Hartree-Fock approximation.In all panels, the experimentally estimated model parameters from [16] are used.Two-magnon scattering dominates the response, producing a streak-like continuum despite the absence of fractionalisation on the α-chains.It is observed that the lateral extent of the rephasing signature is greatly diminished in high fields, due to the energy cost of magnon pair creation increasingly dominating over the width of the 2-magnon continuum.For all figures, L = 100, and the pulse ratio is A0 : Aτ = 1 : 5.
measurement specifically picks out the lower band edge.This polarization choice is then particularly favorable here as a measure of the proximate criticality on the β chains.
In order to calculate the response of the β chains in this setting, one can perform a basis transformation and rotate the spin basis on every second site by π around the y-axis.This returns gives us back an alternating probe field (as for the [1 10] polarization), but with the signs of J x and J z reversed.Looking at Eqs. [( 21)-( 23)], we see that with this change sin ϕ k is very close to vanishing, and is only large near k = π/2, which now corresponds to the bottom band edge.Thus the β chain response is a large spike of intensity close to the origin, which vanishes precisely at J x = J y .
On the α chains, the situation is now somewhat more complicated, as the chains are subject to an external field H which alternates in sign, and a probe field which does not.For a suitable choice of basis, as discussed in Section III A, the excitations on the α chains can be understood to exist in two bands, with the probe field creating pairs of particles of opposite bands, and driving momentum preserving transitions between bands.
In Fig. 12(a), the full 2DCS response for Ce 2 Zr 2 O 7 for a [001] probe field, calculated in the Hartree-Fock approximation, is presented for H = 0.2T and parameters taken from [16].A large central spike of intensity is observed, originating from the β chain contribution to the full response.We find that the β chain contribution is much stronger than that of the α chains, and there is therefore little variation in the full response with field strength.
For small increases in J x − J y , as seen in Fig. 12(b) and (c), the central peak is resolved into the six distinct signatures produced by χ (3) (t, t, t + τ ).This is reflects both the widening of the gap to excitations, which moves each signature away from the origin, and sin ϕ k becoming large for values of k further from π/2.
Approaching the critical point, these six peaks coalesce together, before disappearing entirely precisely at the isotropic point.
While the β chain response vanishes at the isotropic point, the α chain response will remain.We show the α chain response separately for varying values of the anisotropy in Fig. 13.Here, the most notable features lie close to the ω t = 0 axis.Referring to Eq. ( 38), this signature is produced by the term A (2) , whose time dependence is sin(2 λk τ ) cos(∆λ k t), where ∆λ k is the splitting between magnon bands, in the two-band picture described in Section III A. The non-zero ∆λ k shifts the response off of the ω t = 0 axis, producing the slightly curved responses observed in Fig. 13.This curvature could then be used to infer the energy changes in these inter-band transition processes, and thereby deduce the asymmetry in the magnon dispersion about k = π/2.

B. Nd2Zr2O7
Calculating the two dimensional coherent spectroscopic response of Nd 2 Zr 2 O 7 within the CHAIN x phase is made more complicated by the presence of a non-zero mixing angle θ.For the β chains, we can still make calculations within the Hartree-Fock approximation, and use the matrix Pfaffian construction for the expectation values involving S x operators now contributing to the response.For the α chains, we make predictions of the 2DCS response only in the high field limit, wherein we treat excitations about the field polarized state perturbatively.
For the β chains, the Hartree-Fock re-scaled parameters are given by: We observe that at the reported value of θ = 0.98, the rephasing response is seemingly very weak, relative to the rest of the signal.Similar to the situation in Ce 2 Zr 2 O 7 , we attribute this to the strengthening of the ω t = 0 and ω t = ω τ parts of the response for non-zero θ.As seen in Fig. 14 (a) and (b), the rephasing contribution is similar in strength to the other parts of the signal for lower values of the field mixing angle.
As with Ce 2 Zr 2 O 7 , it is only the strengths of different parts of the response that vary with θ.The locations of response peaks are unchanged, which is to be expected, as the β chain Hamiltonian is independent of θ, and it is only the matrix elements with the probe field that are changing.
It is notable that values θ not close to a multiple of π/2 seem to disfavour the rephasing part of the response, as this the most useful signature which cleanly separates out the contributions from fractionalisation and finite lifetimes to the breadth of a continuum [5].Thus whilst the general expectation is that 2DCS enables unambiguous measurements of spinon continua, the finite value of θ for Nd 2 Zr 2 O 7 seems to obscure this information somewhat.This is different to the case of Ce 2 Zr 2 O 7 , discussed above, where θ is small or vanishing and the rephasing signature is predicted to be clear (see Fig. 9).

[110] polarization: 1-and 2-magnon excitations
If we instead probe the α chains of Nd 2 Zr 2 O 7 using a [110] field, we introduce for finite θ a coupling to the S x pseudo spins, or equivalently a non-zero value of J xz , which allows us to probe singly created magnon excitations.As discussed in Section III C, the S x terms in the Hamiltonian become highly non-local in Jordan-Wigner fermions, and so the leading order behaviour of the response is accessed in perturbation theory.The predicted response in the high field limit is given by Eq. (60) and Eq.(61).In general, we predict that both one-and twomagnon excitations will be observable in this polarization, with their relative strength in the signal being controlled by θ.For the value of θ found in [26] for Nd 2 Zr 2 O 7 , the one-magnon contribution is found to dominate.
In Fig. 15, the α chain response in the perturbative limit is plotted for varying values of θ.Fig. 15(a) shows the response for the experimentally estimated parameters, wherein four response peaks are observed.In the high field limit of the exactly soluble model, the rephasing response decays as (J µ /h) 4 , whilst other signatures fall away more slowly with (J µ /h) 2 .We similarly see here that, keeping terms only up to the second power in J µ /h, no rephasing-like signature is observed.
The four peaks that can be seen occur at energies roughly equal to h, the energy cost of a single spin flip, indicating that these are indeed single magnon excitations.As θ is tuned away from its reported value to 0.5 in the middle panel of Fig. 15, the amplitude of the two magnon signature increases, and both sets of peaks can be seen simultaneously, and are distinct from one another.Tuning θ to zero completely, we observe that the single magnon peaks vanish, leaving only the higher energy two magnon signatures.
The response from the two magnon sector is not visible at θ = 0.98 because the value of J xx − J yy is much smaller than J xz for the reported model parameters, and so the probe field has a much greater amplitude to excite magnons one at a time.J xx − J yy , and hence the two magnon response, vanishes completely at θ = ±1.07(2)and θ = ±2.07(0)radians.Tuning θ towards 0, J xx − J yy increases whilst J xz decreases with sin(2θ).This leads  to the observed cross-over between a dominant single magnon response and a dominant two magnon response.

3.
[001] polarization: observing rephasing at large θ Subjecting Nd 2 Zr 2 O 7 to a [001] probe produces responses from both α and β chains.The contribution to the response from the β chains is qualitatively very similar as for the [1 10] polarization, whilst the α chain response shows signatures of probe field induced transitions between k = 0 and k = π states.The α chains are again analysed in a high field perturbative limit, whilst the β chain response is calculated using the matrix Pfaffian expansion in the Hartree-Fock approximation.As differing methods are used for each type of chain, the responses for each are presented separately.
First, as discussed in Section III C, we find that the α chain response comes only from single magnon excitations to second order in perturbation theory.Fig. 16 presents the predicted response for values of θ varied around that given in [26].Response peaks occur at the excitation energies of k = 0 and k = π magnons, as well as at the transition energy ∆λ = λ 0 − λ π .As in Section IV A 3, this can be understood in a two band picture with the doubling of the unit cell coming from the different position dependencies of the static external magnetic field and the probe fields.It is observed that the two peaks on the line ω t = 0 each split in two, and move apart perpendicular to this axis a distance ∆λ.Also, a further set of peaks of opposite sign appear near the origin, which vanish if ∆λ = 0.
As the mixing angle θ is varied, the difference between the k = 0 and k = π magnon energies changes also, and   [26].We observe that θ controls the relative visibility of the one-and two-magnon contributions to the response.For θ = 0.5 both contributions are observable, leading to two distinct sets of peaks.The lower energy, single magnon, peaks are also visible at θ = 0.98.At θ = 0, only the two magnon response is visible.For all figures, L = 40, and the pulse ratio is A0 : Aτ = 1 : 5.
we find that as the mixing angle is increased from π/4, through θ = 0.98, to θ = 1.2, the increasing value of ∆λ can be observed in the increasing separation of peaks either side of the ω t = 0 axis.The contribution from the β chains is presented in Fig. 17, as calculated using the matrix Pfaffian construction.The response is calculated for various values of θ, tuning the matrix elements between the pseudospins and the probe field.The most noticeable difference seen in the [001] β chain response when compared the [1 10] response, is that at the reported mixing angle of θ = 0.98 [26] the rephasing signal is far more visible, and is not completely dominated by other signatures in the response.Examining Fig. 17, one can see that the relative intensities do still evolve with θ for this polarization choice, with the relative amplitude of the rephasing signature being higher in Fig. 17    than at θ = 0.98.We note that the total intensity of the response from the β chains using a [001] probe is somewhat weaker than for the [1 10] polarization.
These results again show that, whilst a large non-zero θ may obscure the rephasing response for the β chains, it is not absent entirely.Further, as it is only the overall intensity that is modulated with θ, this response can still be used to probe properties of fractionalized spinon excitations on the β chains.They also suggest the [001] polarization as more favourable than [1 10] for the excitations of the β chains in Nd 2 Zr 2 O 7 .Plots are generated in the Hartree Fock approximation, making use of the matrix Pfaffian expansion, which introduces a global rescaling of the response due to approximations made in computations.Much as in Fig. 14, only variation in the relative intensities of different parts of the response is observed as θ is changed.Notably, the rephasing signal is not so heavily suppressed at θ = 0.98 as it was for the [1 10] probe field, as seen in (b), though the ωt = 0 signal again dominates if θ is further increased, as seen in (c).In all cases, L = 40, and A0 : Aτ = 1 : 5.

V. CONCLUSION
In this article we have explored how the polarization of magnetic probe fields in two dimensional coherent spectroscopy influences the measured response, and we have demonstrated how this polarization can be used as a powerful control parameter.By analysing the detailed coupling as a function of polarisation, we can infer information ranging from the detailed nature of the microscopic degrees of freedom all the way to the determination, and isolation, of different types of excitations in an exotic magnet.
We have focused on the materials Ce 2 Zr 2 O 7 and Nd 2 Zr 2 O 7 , both fascinating examples of frustrated quantum magnets presenting field induced Chain phases of effectively one dimensional character [21,22].By vary-ing the probe field polarization, the 2DCS response from fractionalised, unpolarised β chains can be isolated from the response of non-fractionalised, polarised α chains.
The strong dependence of the 2DCS response on the coupling between the probe field and a system's degrees of freedom is also highlighted in how the dipolar-octupolar mixing angle θ influences the response, with this mixing angle being observed to control the relative intensities of one-and two-magnon signatures.
In the particular case of Ce 2 Zr 2 O 7 , θ = 0 and we observe that the response to a [1 10] probe, which accesses the β chains, is qualitatively similar to the non-fractionalised response of the α chains to a [110] field, as the magnetic dipole operator is only able to create magnons in pairs.In Nd 2 Zr 2 O 7 , θ is non-zero, and single magnon excitations dominate the [110] response.We see also that large values of θ can somewhat suppress the relative amplitude of the rephasing streak along the ω t = −ω τ line.
A different polarisation choice can be used to probe another aspect of the physics of Ce 2 Zr 2 O 7 .The chain phase of this material is believed to be close to a critical point, and we find that the [001] polarization serves as a particularly sensitive measure of the proximity to criticality.
We hope that the predictions presented here will inspire experiments to test them, exploring the quantum excitations of field induced chain phases.This may be challenging on a technical level, since the relatively weak exchange interactions in these systems (∼ 0.1meV) mean that both excellent frequency resolution and very low temperatures would be required.Nevertheless, there is no reason of principle as to why these obstacles cannot be overcome and we expect that as 2DCS becomes a more widely used tool in the study of quantum magnetism, the approach will be increasingly optimised to meet these challenges.

FIG. 1 .
FIG. 1.(a) The pyrochlore lattice and its decomposition into α and β chains in a [110] magnetic field.α chains run parallel to the field, and are coloured red in the left panel.In a sufficiently large applied field, spins on the α chain become polarized, whilst the perpendicular β chains do not feel the field.(b) A single tetrahedron demonstrating one possible alignment of spins that satisfy the '2-in-2-out' rule while also polarising the α chain spins with respect to the applied field.The β spin chains can then be in one of two configurations that preserve the ice-rule.The numerical labels 0-3 indicate the labelling convention for the four fcc sublattices Lµ used throughout this Article.
Phases in Dipolar-Octupolar Pyrochlores in a 110 magnetic field

FIG. 2 .
FIG. 2. Ground state phase diagrams for Ce2Zr2O7 for a finite [110] field, for parameters in the vicinity of the values determined in [16].Phase diagrams are shown as a function of field strength and γ = Jx − Jy.In all figures, the two vertical dotted red lines represent the predicted model parameters of Ce2Zr2O7 with increasing magnetic field H, which place the material close to a transition between states with dominant correlations along the x and ỹ pseudospin axes.(a) Classical mean field phase diagram with θ = 0 [Eq.(3)-(5)].In this case, with θ precisely vanishing, the system transitions to a CHAIN * λ regime at finite field (λ = x, ỹ), and then through a second order transition to a CHAIN λ phase for a critical value of H. (b) Classical mean field phase diagram with θ = 0.01.Here the transition at finite field vanishes for Jx > Jy(γ > 0), due to the small but finite value of θ [30].(c) Phase diagram including the effects of intra-chain quantum fluctuations with θ = 0.It is observed that the CHAIN λ phases, within which the α chains are polarized parallel to the [110] field, are greatly stabilised, with the CHAIN * λ Let us specialise now to Ce 2 Zr 2 O 7 and Nd 2 Zr 2 O 7 .These materials have the following experimentally estimated model parameters: (J x , J y ) = (0.064, 0.063) meV J z = 0.011 meV; θ ≈ 0 rad; g z = 2.57, for Ce 2 Zr 2 O 7 [16, 21], with (...) indicating that the values of J x and J y can be swapped without affecting the agreement with experiment , and J x = 0.091(9) meV; J y = 0.014(6) meV; J z = −0.046(2)meV; θ =0.98(3) rad; g z = 5.0(1), compare to χ (n) µ,direct greatly alters the resulting 2DCS response, a fact that makes the [001] polarization channel particularly useful in analysing Ce 2 Zr 2 O 7 and the proximity of its model parameters to the critical point.B. Effect of Jz ̸ = 0: Hartree-Fock Approximation Thus far we have only calculated the required non-linear magnetic susceptibilities, χ (n) µ,alt and χ (n) (n) β,alt and χ (n) β,direct , the effects of finite

FIG. 4 .
FIG. 4. Parameter flows under the Hartree-Fock approximation as Jz is increased for the zero field case ((a)), and for an initial field h = 1 ((b)).The blue dots show starting parameter values, and the curves show the evolution of the re-scaled parameters as Jz is increased from 0 to 1.The red dashed lines indicate phase boundaries.The central region of (b) is the field polarized paramagnetic regime, the four other phases are ferro-or antiferromagnetic phases with dominant correlations along either x or ỹ spin axes.In (b), note that the axes are re-scaled to Jx/h and Jy/h, such that the phase diagram can be presented in two dimensions.We observe that no phase boundary crossings occur with increasing Jz in the zero field case, whilst finite Jz can generate transitions into the central paramagnetic phase in the presence of a finite field.

FIG. 6 .
FIG.6.Evolution of re-scaled Ce2Zr2O7 model parameters in the Hartree-Fock approximation with varying strength of the applied field H.We observe only small modifications to all parameters, with H ′ scaling approximately linearly with H.It is also noted that the anisotropy Jx−Jy remains approximately constant with increasing H.

FIG. 7 .
FIG. 7. Demonstration of the robustness of the rephasing signal ((b) and (d)) to broadening effects, when compared to the ωt = ωτ signature ((a) and (c)).(a) A set of nine peaks along the ωt = ωτ with small broadening Γ = 0.1.The function varies rapidly from positive to negative between each peak.(b) A similar response from nine closely spaced peaks, but this time organised along the ωt = −ωτ line.Here the function is always negative along the length of the response, and is flanked by positive regions.When broadening is increased to Γ = 1, we see in (c) that a broad signal is obtained, whilst in (d) a diagonal streak is still clearly observed along the ωt = −ωτ line.

FIG. 9 .
FIG. 9. 2DCS response of β chains in Ce2Zr2O7 to a [1 10] polarized probe field, for varying values of the anisotropy Jx − Jỹ.(a) Jx − Jỹ = 0.001meV, as in the best fit parameters from[16], (b) Jx − Jỹ = 0.014meV, (c) Jx − Jỹ = 0.027meV.Calculations are made using Eq.(25), within the Hartree-Fock approximation on a system of 100 spins, along with large phenomenological depopulation and dephasing times T1 = T2 = 150meV −1 .The energy gap is observed to grow as the β chain Hamiltonian is taken further away from criticality, with the low energy edge of the rephasing streak along ωt = −ωτ moving away from the origin.The intensity of the rephasing streak closely tracks the spinon density of states, as shown in Fig.8.The strength of the magnetisation response is in arbitrary units, and the ratio of the strength of the first field pulse to the second is 1 : 5.

3 .
[001] polarization: probing nearby criticality Finally, we consider what the 2DCS response using a [001] polarized probe field can tell us about the excitations of Ce 2 Zr 2 O 7 .Such a probe field couples to both α and β chains, and does so uniformly across each chain, instead of alternating in sign from one site to the next.The relative π phase shift in the momentum of the probe field to the [110] and [1 10] probe fields (and to the external field H) leads to qualitatively different 2DCS responses.The response of the β chains to the [001] polarization is very different to the case of the [1 10] polarization: instead of being sensitive to the full spinon density of states, the

FIG. 10 .
FIG.10.2DCS response of β chains in Ce2Zr2O7 to a [1 10] polarized probe field, for increasing values of the dipolar-octupolar mixing angle θ.Responses are generated using the matrix pfaffian method, and so the overall scale of plots in this figure differ from those generated without this expansion.(a) θ = 0, as in the best fits in[16], (b) θ = 0.3, (c) θ = 0.5.We observe that as θ is increased from zero, a large signal is introduced along the ωt = 0 line, which dominates over the rephasing signal.A chain length of L = 40 was used for these calculations, with A0 : Aτ = 1 : 5.

FIG. 11 .
FIG. 11. 2DCS response of α chains in Ce2Zr2O7 to a [110] probe field, for increasing strength of the external field H. Figures are produced using Eq.(25), within the Hartree-Fock approximation.In all panels, the experimentally estimated model parameters from[16] are used.Two-magnon scattering dominates the response, producing a streak-like continuum despite the absence of fractionalisation on the α-chains.It is observed that the lateral extent of the rephasing signature is greatly diminished in high fields, due to the energy cost of magnon pair creation increasingly dominating over the width of the 2-magnon continuum.For all figures, L = 100, and the pulse ratio is A0 : Aτ = 1 : 5.
. [1 10] polarization: absence of spinon-rephasing signal for large θ In Fig. 14, we present several figures showing the 2DCS response of the β chains in Nd 2 Zr 2 O 7 to a probe field in the [1 10] direction, for the experimentally estimated model parameters, with varying θ.

FIG. 12 .
FIG.12.Total 2DCS response from α and β chains in Ce2Zr2O7 for a probe field in the [001] direction, calculated in the Hartree-Fock approximation, while varying Jx − Jỹ.A field strength of H = 0.2T is used in all panels.(a) Jx − Jỹ = 0.001meV as in the best fits from[16], (b) Jx − Jỹ = 0.014meV, (c) Jx − Jỹ = 0.027meV.The central peak is observed to resolve into the six distinct signatures produced by χ (3) (t, t, t + τ ), which dominates for the pulse strength ratio chosen of A0 : Aτ = 1 : 5.As the anisotropy is increased, the β chains are driven further from criticality, and the six features spread out from the center of the plot.At the isotropic point, the β chain response collapses completely to ωt = ωτ = 0, and vanishes, leaving only the much weaker α response.

FIG. 13 .
FIG.13.α chain contribution to the 2DCS response of Ce2Zr2O7 to a [001] field, calculated in the Hartree-Fock approximation, while varying Jx − Jỹ.A field strength of H = 0.2T is used in all panels.(a) Jx − Jỹ = 0.001meV as in the best fits from[16].(b) Jx − Jỹ = 0.014meV.(c) Jx − Jỹ = 0.027meV.We observe that the α chain response is roughly an order of magnitude weaker than the β response.As the anisotropy is increased from left to right, we observe that the curvature in the signature close to the ωt = 0 axis increases, indicating the contribution of inter-band magnon transitions driven by the probe field.

FIG. 14 .
FIG. 14. 2DCS response of β chains in Nd2Zr2O7 to a [1 10] field, for varying values of the parameter θ.The other parameters take their reported values from[26].(a) θ = 0, (b) θ = 0.5 and (c) θ = 0.98 as found in[26].Calculations are made using the Hartree-Fock approximation, and by making use of the matrix Pfaffian method, which introduces an order one global rescaling of the responses due to approximations used in computations.Both the overall intensity of the response, and the relative intensity of different signatures vary with θ, with a noticeable shift in relative intensity away from the rephasing signal, but no other qualitative changes to the response are observed for the β chains over this range of dipolar-octupolar mixing angles.In all cases, L = 40, and A0 : Aτ = 1 : 5.

FIG. 15 .
FIG.15.2DCS response of α chains in Nd2Zr2O7 to a[110]  probe field, as calculated using perturbation theory in the high field (large h) limit, with varying mixing angle θ.Parameters apart from θ are set to their values from[26].(a)θ = 0.0, (b) θ = 0.5, (c) θ = 0.98 as in[26].We observe that θ controls the relative visibility of the one-and two-magnon contributions to the response.For θ = 0.5 both contributions are observable, leading to two distinct sets of peaks.The lower energy, single magnon, peaks are also visible at θ = 0.98.At θ = 0, only the two magnon response is visible.For all figures, L = 40, and the pulse ratio is A0 : Aτ = 1 : 5.

FIG. 16 .
FIG.16.2DCS response of only the α chains in Nd2Zr2O7 to a [001] probe field, as calculated using perturbation theory in the high field limit, with varying mixing angle θ.Parameters apart from θ are set to their values from[26].(a)θ = π/4, (b) θ = 0.98 as in[26] and (c) θ = 1.2.We observe that as θ is increased over this range of values, the two peaks on the ωt = 0 line separate into two signals, and a third set of peaks is generated close to the origin.This splitting is due to transitions between k = 0 and k = π magnon states induced by the probe field, which has an effective momentum of π in the basis in which the field polarized ground state is uniformly ordered.

FIG. 17 .
FIG. 17. 2DCS response of just the β chains in Nd2Zr2O7 to a [001] field, for varying values of the parameter θ.Parameters other than θ take their reported values[26].In (a) θ = π/4, in (b) θ = 0.98, the best fit value from[26], and in (c) θ = 1.2.Plots are generated in the Hartree Fock approximation, making use of the matrix Pfaffian expansion, which introduces a global rescaling of the response due to approximations made in computations.Much as in Fig.14, only variation in the relative intensities of different parts of the response is observed as θ is changed.Notably, the rephasing signal is not so heavily suppressed at θ = 0.98 as it was for the [1 10] probe field, as seen in (b), though the ωt = 0 signal again dominates if θ is further increased, as seen in (c).In all cases, L = 40, and A0 : Aτ = 1 : 5.