Thermal masses and trapped-ion quantum spin models: a self-consistent approach to Yukawa-type interactions in the λϕ 4 model

,


Introduction
The field of quantum technologies, which aims at developing quantum devices that provide novel functionalities with a quantifiable advantage with respect their classical counterparts, has become a promising area of research in both the academic sector and the technological industry, e.g.[1,2].Regarding the application of these technologies to quantum computation [3], there has been a remarkable recent progress [4][5][6][7][8] towards the long-term goal of a large-scale fault-tolerant error-corrected device [9] that can outperform classical computers in relevant tasks [10].However, before these large-scale devices become available, one is restricted to operate with small to mid-scale prototypes in the so-called noisy intermediate-scale quantum (NISQ) era [11].Here, one aims at realizing specific circuits, or even prototype quantum algorithms [12], on the largest possible number of qubits and gates, evading the overhead of active quantum error correction.Remarkably, even in the presence of noise, these NISQ devices have already enabled the demonstration of quantum advantage [13,14], that is, using a quantum device to solve a problem that would require an unfeasible amount of time with a classical machine.A current research goal is to extend these demonstrations of quantum advantage to problems that can be of practical relevance in various areas of science.
In this respect, the simulation of quantum manybody models is a problem of considerable interest in different disciplines, ranging from quantum chemistry, to condensed matter and high-energy physics.As originally emphasised by Richard Feynman [15], the characteristic exponential scaling of the size of the Hilbert space of a quantum many-body system hints at the inherent complexity of this type of problems, and the inefficiency of a brute-force numerical approach based on classical computers.Although various numerical methods have been developed over the years to overcome these difficulties, there are still many open questions regarding real-time dynamics, finite-fermion densities and, generally, stronglycorrelated phenomena.The idea of quantum simulations [16][17][18] is to use a quantum device instead of a classical one, which can be controlled so as to reproduce the properties and dynamics of the model of interest.This has already found several applications in the aforementioned areas [19][20][21][22][23].
A quantum simulation proceeds by, first, encoding the degrees of freedom of the target model into those of the quantum device, preparing a specific initial state, and then letting the system evolve in time before a measurement stage.This can be achieved in two ways.One can use the same building blocks as in quantum computers, i.e. qubits, which are then acted upon by a sequence of quantum logic gates.This sequence reproduces approximately the real-time dynamics of the model under a Suzuki-Trotter expansion [24], and leads to the so-called digital quantum simulations [25].Note that, in spite of working with qubits, one can simulate fermionic and bosonic degrees of freedom at the expense of an overhead in the number of gates and/or qubit-register size, e.g.[26].Alternatively, one can use special-purpose quantum simulators that already have spins, fermions, bosons, or combinations thereof, as the relevant degrees of freedom.This advantage comes at the price of certain limitations in the range of models that can be simulated since, in general, one cannot realize an arbitrary unitary on the exponentially-large Hilbert space.In fact, these quantum simulators are not acted upon by concatenating gates drawn from a universal gate set, but rather by letting the system evolve continuously in time under approximate effective Hamiltonians with a restricted set of terms.By tuning the strength of these terms, one can mimic approximately the target model in a specific parameter regime.These devices are known as analog quantum simulators [27].
In this manuscript, we are interested in the use of long trapped-ion chains as analog quantum simulators [21,22], targeting models from high-energy physics.The potential advantage of working with a large number of laser-cooled ions for frequency standards and clocks led to the development of linear Paul traps [28], which can store long ion chains along the trap symmetry axis while reducing the micromotion.In the digital approach, these chains are operated by laser/microwave radiation to implement a highfidelity gate set, which has been exploited for quantum simulations of small-scale spin models in condensed matter, either following a Hamiltonian [29,30] or a Lindbladian [31,32] time evolution.Additionally, digital quantum simulations of a lattice gauge theory [33][34][35] have also been performed [36,37].
For further scaling, we consider in this work the analog quantum simulation approach.Here, instead of using sequences of single and two-qubit gates [38,39], one can obtain an effective spin model for the whole ion chain by acting with always-on far-detuned lasers, which typically lead to long-range spin-spin interactions mediated by the phonons of the chain [40].This idea has turned out to be extremely fruitful for the analog quantum simulation of magnetism [41][42][43][44][44][45][46][47][48][49][50][51].In all of these experiments, the focus lies on the interacting spins, while the phonons are mere auxiliary degrees of freedom.In contrast, as advocated in recent works [52,53], including these bosonic degrees of freedom in the quantum simulation offers a more complete picture, and provides a neat understanding of some characteristics of the spin models.In particular, the decay of the spin-spin couplings with the inter-ion distance can be described by a dipolar decay with an additional exponential tail that is controlled by the value of the laser detuning [54,55].The origin of this particular distance dependence is clarified by a long-wavelength description of the model [53], in which the phonons are described as quantised sound waves in terms of a relativistic Klein-Gordon field ϕ(x) [56,57].This scalar field has a Yukawa-type coupling to the spins, and can be used to mediate the spin-spin interactions and explore the entanglement dynamics [53,58].
An interesting question from the perspective of QFTs is to push this analogy further, and move to situations beyond the free Klein-Gordon QFT.Inspired by the fermion-Higgs sector of the electroweak theory [59], our trapped-ion quantum simulation would become richer in the presence of λϕ 4 interactions [60].As discussed in [53], this type of self-interactions be-come the relevant ones in the vicinity of a structural phase transition [61], and lead to a variety of scattering events for the bosonic excitations of the scalar field.These quantum effects can be expressed in terms of Feynman loop diagrams leading to changes in the physical mass of the bosons.Interestingly, in the trapped-ion context, this renormalization effect can affect the effective spin-spin interactions of the quantum simulator as one approaches the structural phase transition, opening an original route for probing the underlying effective QFT via the real-time dynamics of the spins.Unfortunately, for typical realizations, the rigidity of the ion chain tends to mask these quantum effects [53,62], such that the change of the physical mass with the quartic coupling λ is predicted to be very small.In this work, we argue that one can overcome this limitation by working at finite temperatures which, in QFTS, lead to the so-called thermal masses and the phenomenon of symmetry restoration [63,64].We present a non-perturbative self-consistent approach to predict the renormalization of the spin-spin interactions at finite temperatures, and discuss how these can be used to probe the thermal QFT that controls the structural phase transition of the ion crystal at long wavelengths.Our presentation is organised as follows.
In Sec. 2, we discuss how to describe a structural phase transition in a trapped-ion chain by an effective λϕ 4 QFT, and how the trapped-ion quantum simulators of spin models can be understood in light of this QFT as a Yukawa-type problem of interactions mediated by a scalar field.We finish by discussing qualitatively the effects of self-interactions of the scalar field and non-zero temperatures in these Yukawa-type spin-spin interactions.In Sec. 3, we derive a set of equations to deal quantitatively with self-interactions and non-zero temperatures in this field theory, which requires resuming certain types of Feynman diagrams in a non-perturbative approach.We emphasise how to deal with ultraviolet and infrared divergences in this approach, both of which arise for the low spacetime dimensions relevant to the trapped-ion problem.In Sec. 4, we discuss in detail how to adapt these techniques to the specific details of trapped ions, and describe our numerical approach to solve the aforementioned equations.Finally, we discuss how these results predict a temperature-dependent contribution to the physical mass of the scalar particles, which changes the range of the Yukawa-type spin-spin interactions, giving concrete predictions for a realistic trapped-ion experiment.Finally, in Sec. 5, we present our conclusions and outlook. 2

Quantum simulation of Yukawa-type models
For the shake of completeness, and to fix our notation, we review in this section some concepts about trapped ions, as well as the connection of the phonon-mediated spin-spin interactions in the vicinity of a structural transition to the quantum simulation of a Yukawatype QFT at non-zero temperatures.

From trapped ions to λϕ 4 quantum fields
As advanced in the introduction, we consider a system of N 1 atomic ions of charge e confined in a linear Paul trap [65,66], and assume that the symmetry axis of this trap lies along the x direction.In this kind of devices, ions are trapped using a combination of AC and DC potentials, and can reach a stable crystalline distribution for hours or even days, depending on the trap design, and the cooling and vacuum conditions [67].In the pseudo-potential approximation [28,68], the secular motion of the ions can be described by an effective quadratic potential with constant trap frequency {ω α } that aims at confining the ions along each axis α ∈ {x, y, z}.Since the ions are charged, there is a competition between this overall trapping potential and the inter-particle Coulomb repulsion, leading to ) where r i , p i are the canonical position-momentum operators of the ions with mass m a and charge e, and ϵ 0 is the vacuum permittivity.This competition leads to a set of equilibrium positions {r 0 i = l r0 i } N i=1 , where we have introduced a constant with the units of length l that fulfills l 3 = e 2 /4πϵ 0 m a ω 2 x [69,70], which are obtained by solving the non-linear equations where κ α = (ω x /ω α ) 2 .For ω x ≪ ω y , ω z , these equilibrium positions lie along the symmetry axis r 0 i = x i e x , and present a typical lattice spacing that is rather homogeneous in the bulk of the ion chain After expanding the Coulomb potential in Eq. (1) to quadratic order in the small displacements r i (t) = r 0 i + u i (t), the secular ion motion is accurately described by H ≈ H 2 with (3) Here, π i,α = m a ∂ t u i,α are the momenta conjugate to the displacement operators [u i,α (t), π j,β (t)] = iℏδ i,j δ α,β , and we have introduced the effective spring constants According to this approximation, the small displacements of the ion crystal are governed by a model that resembles a simple elastic/harmonic chain [71].In comparison to this textbook example, we note that the effective spring constants do not restrict to nearest neighbors and, moreover, one also finds an additional local elastic term k α = m a ω 2 α .For the longitudinal displacements u i,x , these differences do not introduce qualitative changes in the physics: the ions vibrate in the form of collective excitations in analogy to the acoustic phonons in a solid, which are associated to quantised compressional sound waves.In contrast, the different sign of the spring constants (4) that couple to the transverse displacements u i,z can actually lead to a situation that differs markedly from the quantised shear sound waves of an elastic solid.In fact, as one rises the ratio of the Paul-trap frequencies κ z = (ω x /ω z ) 2 above a critical value κ z > κ z,c (N 1 ) [72,73], there is a structural transition from the linear chain into what is known as a zigzag ladder, as first observed experimentally in [28].For a specific ion number and frequency ratio, one can numerically solve the system of non-linear equations for the equilibrium positions (2), leading to the configurations displayed in Fig. 1.The critical value decreases with the number of ions as a power law κ z,c (N 1 ) = aN b 1 for certain constants a > 0 and b < 0 [72]; a scaling verified in experiments [74].
These structural changes are the finite-size precursors of a quantum phase transition, which can be characterised by the spontaneous breaking of a discrete inversion symmetry with respect to the trap axis.Even working far away from the thermodynamic limit N 1 → ∞, these mesoscopic structural transitions have a well-defined soft mode [61], which allows one to derive a low-energy/long-wavelength approximation that goes beyond the elastic limit.This requires going beyond the quadratic Hamiltonian in Eq. (3) by considering higher orders in the expansion of the Coulomb interaction (1).This long-wavelength theory is the counterpart of the λϕ 4 model QFTs in which the aforementioned inversion symmetry is ϕ(x) → −ϕ(x) and gets spontaneously broken at a certain critical point.In the Hamiltonian formulation, this QFT can be written as (5) where x = (t, x) are the Minkowski spacetime coordinates, and m 0 , λ 0 are the bare mass and bare coupling constants, respectively.In this expression, we have used natural units ℏ = c = 1, as customary in high-energy physics.For the connection to trapped ions, it is more appropriate to work in SI units [75], such that (6) In comparison to natural units, where the field operator and its conjugate momentum have scaling dimensions [ϕ] = L 0 , [π] = L −1 , the dimensional analysis in SI units leads to the following scaling Let us discuss first how the quadratic Klein-Gordon part of the QFT (6) arises through a coarse-graining procedure of Eq. (1), when focusing on the quadratic approximation (3).In this way, we will highlight the important differences with respect to the effective QFT for compressional sound waves in an elastic chain [71].By focusing on the bulk of the ion crystal, we can use periodic boundary conditions, and move to a Fourier representation for the small displacements Here, the quasi-momentum is k = 2π dN1 n 1 for n 1 ∈ {1, • • • , N 1 }, and thus lies within the first Brillouin zone BZ = (0, 2π/d], where d is the lattice spacing in the bulk of the ion chain, which is approximately constant (see Fig. 1 (a)).Using this transformation, and focusing first on the quadratic part (3), one can derive the dispersion relation This dispersion relation resembles that of a standard lattice regularization of the Klein-Gordon QFT in Eq. (6).In a Hamiltonian lattice field theory [34], the spatial components are discretised as x → x i = ia, where i ∈ {1, ..., N 1 }, which require introducing an artificial lattice spacing a that regularises the QFT by an ultra-violet cutoff Λ c = π/a.Following the lattice approach [76], one performs the following substitutions (9) By applying a similar Fourier transformation (7) to the discretised QFT (6), one finds the dispersion relation where Although this expression clearly resembles the trapped-ion case (8), there are important differences.The most apparent one is that the dispersion relation in Eq. (8) contains a dipolar tail due to the long-range nature of the effective spring constants (4).In addition, there is a sign difference with respect to (10) that will be crucial for the long-wavelength coarse-graining.
The dispersion relation of the Klein-Gordon QFT ω 2 (k) = m 2 0 c 4 /ℏ 2 + c 2 k 2 is recovered from Eq. (10) by considering |k| ≪ Λ c .This particular low-energy limit corresponds to a coarse-graining approximation, where the relevant scale is much larger than the lattice spacing ξ 0 ∝ 1/m 0 ≫ a, and one says that the continuum QFT is recovered by sending a → 0. The coarse graining in the trapped-ion case (8) is slightly different.Due to the sign difference in Eq. (8), the lowest-energy mode corresponds to the so-called zigzag mode, and one has to expand around k = π/d + δk.From this expansion, one can readily recover the same Klein-Gordon dispersion, identifying the speed of the transversal sound waves as which plays the role of an effective speed of light in the quantum simulator of the QFT (6).In this expression, we have used a truncated version of the Dirichlet eta function, namely such that η N1 (1) → log 2 in the thermodynamic limit N 1 → ∞.In addition, one obtains the effective bare mass m 2 0 =: where ω zz is the frequency of the zigzag mode, and we have introduced a truncated version of the Riemann zeta function such that ζ N1 (3) → 1.202 in the thermodynamic limit N 1 → ∞, corresponding to Apéry's constant.At this point, we emphasise that the coarse-graining procedure in a physical trapped-ion lattice of spacing d does not require sending the spacing to d → 0. Alternatively, one can tune the parameters close to a critical point where the bare mass (13) would vanish m 0 → 0, such that the effective Compton wavelength fulfills ξ 0 = ℏ/m 0 c t ≫ d, and the low-energy physics does not depend on microscopic lattice details, but rather on the universal properties captured by a QFT.As it can be checked from Eq. (13), this critical point m 2 0 | c = 0 coincides precisely with the linear-to-zigzag transition at This critical point κ z,c (N 1 ) has a scaling with the number of ions that agrees with the previouslymentioned power laws [72,73], already for moderate ion numbers [77].So far, our discussion has focused on the elastic/quadratic terms and the dispersion relation, but we have not identified the fields yet.In order to find their trapped-ion analogue, we need to separate the rapid oscillations around the zigzag mode [53,78] from a slowly-varying envelope that will play the role of the coarse-grained field [71,79].In addition, one has to rescale the position and momentum operators in Eq. (3) to achieve the correct scaling dimensions below Eq. ( 6) ) such that one recovers the canonical algebra [ϕ(t, x), π(t, y)] = iℏδ i,j /d → iℏδ(x − y).Since the coarse-grained fields vary slowly, one can perform a gradient expansion ϕ(t, y) ≈ ϕ(t, x) + (y − x)∂ x ϕ(t, x), and obtain the Klein-Gordon part of Eq. (6) from the original microscopic theory (3).
Let us note that, at the level of the effective QFT (6), the critical point m 2 0 = 0 is stable thanks to the additional quartic potential with λ 0 > 0. In the trapped-ion case, one needs to extend the microscopic theory (3) to quartic order [70], and apply again the gradient expansion to identify the analogue of the bare quartic coupling.In this procedure, the rigidity of the trapped-ion chain becomes relevant.As occurs with other long-wavelength descriptions in condensedmatter and high-energy physics [80][81][82], one can use thermodynamic arguments when constructing the effective field theory.In the present case, where the coarse-grained QFT can be directly obtained from the trapped-ion model via the gradient expansion, one finds that, in addition to the effective speed of light (11) and bare mass (13), the rigidity modulus of the trapped-ion chain [83,84] gives rise to an addi-tional dimensionless Luttinger paramater [53], namely This parameter quantifies the rigidity of the trappedion chain under a shear strain that aims at deforming it transversely, and becomes important when identifying the coupling constant of the λϕ 4 model.In fact, after expanding Eq. (1) to fourth order and identifying the terms that are more important for the structural phase transition, one finds H ≈ H 2 + H 4 where and we have introduced the quartic coupling matrix At this point, one can perform again the gradient expansion below Eq. ( 16), which allows one to find the final microscopic expression for the λϕ 4 coupling We have thus discussed how a trapped-ion chain in the vicinity of a structural phase transition serves as a quantum simulator of a regularised self-interacting QFT (6) with the effective speed of light in Eq. (11), and the bare parameters in Eqs.(13), (17) and (20).In this context, note that the critical point (15) is obtained by setting the bare mass (13) to zero, and thus corresponds to the classical field-theory calculation in which the minimum of the quartic potential underlying Eq.(6) changes from a single to a double well.At this point, the Z 2 inversion symmetry of the real scalar field ϕ → −ϕ gets spontaneously broken.From the perspective of QFTs, one knows that the classical quartic potential gets quantum contributions in terms of Feynman loop diagrams [85,86], such that the single to double-well will no longer be characterised by the classical critical point.Instead, the excitations are dressed leading to a physical mass m 2 0 → m 2 P , and one finds that the phase transition m 2 P = 0 yields a critical line in parameter space (m 2 0 , λ 0 ) that separates the symmetry-broken and symmetry-preserved regions.Going back to the trapped-ion problem, the critical point (15) will flow with the coupling strength κ z,c → κ z,c (λ 0 ), defining a line that separates the linear chain κ z < κ z,c (λ 0 ) from the zigzag ladder κ z > κ z,c (λ 0 ).
Note that, if one takes the continuum limit a → 0 in the lattice field approach (9), the UV divergences of the loop integrals must be subtracted from the bare mass in order to get finite parameters and draw a meaningful phase diagram.In the trapped-ion case, on the contrary, the lattice spacing remains constant as one approaches the critical point, and it is the physical Compton wavelength which becomes very large ξ 0 → ξ P ≫ d, justifying the long-wavelength description.Accordingly, one can stick to the bare parameters without any additional subtraction, and still find a meaningful phase diagram.One must bear in mind that the critical line will depend on the lattice spacing and other non-universal microscopic properties.In contrast, as one approaches this critical line, the universal properties of the phase transition, i.e. scaling critical exponents, should be controlled by the fixed point of the continuum QFT (6).This corresponds to the so-called Wilson-Fisher fixed point, which can be characterised by an ϵ-expansion in higher dimensions [87].In D = 1 + 1 dimensions, however, the perturbative renormalization-group techniques [60] underlying this ϵ-expansion break down, as all perturbations are relevant in the renormalization-group sense.Localising the critical line of the lattice model, as well as the critical exponents of the corresponding fixed point, requires using non-perturbative techniques, such as Monte Carlo or tensor-network methods [88][89][90][91][92][93][94][95][96].The continuum limit of these studies is consistent with a QFT presenting a second-order phase transition where the Z 2 symmetry gets broken [97], and where the fixed point lies in the universality class of the two-dimensional Ising model.

Yukawa-type interactions and real-time spin dynamics
As advanced in the introduction, we are interested in the use of trapped-ions as quantum simulators of Yukawa-type spin models, and how the real-time dynamics of the spins can serve to probe the underlying interacting QFT.So far, however, we have only discussed the motional degrees of freedom of the ions, leading to the effective QFT (6) in the vicinity of the linear-to-zigzag transition (15).As noted in the introduction, the ions also have an atomic structure with many electronic levels, among which one can select a pair of long-lived states to encode the degrees of freedom of a chain of spins {|↑ i ⟩ , |↓ i ⟩} N1 i=1 .These two states can correspond to the so-called optical, hyperfine or Zeeman qubits in trapped-ion quantum computing [66].Building on this seminal proposal, the subsequent experimental and theoretical efforts have turned trapped ions into one of the leading platforms in the quest of building a large-scale faulttolerant quantum computer [66].In addition to the experiments that we have already mentioned [4,7], which contribute to the continuous effort [98][99][100][101][102][103][104][105][106][107][108][109][110] towards trapped-ion quantum error correction [111], a variety of NISQ algorithms have also been realized over the years [112][113][114][115][116][117][118].The success of these implementations relies on the very accurate performance of the universal gate set [119][120][121][122][123][124].This also includes high-accuracy two-qubit gates which, typically, are the most challenging part of the gate set in any platform.Moreover, either by exploiting ion shuttling [125][126][127] or individually-addressed laser beams [128][129][130], one can implement these entangling gates between arbitrary qubit pairs in the crystal.This allows for a programmable connectivity that has allowed to demonstrate complex protocols with current NISQ error rates, e.g.[4,7], which would not be possible if the connectivity was local.We note that the digital methods can also be combined with classical variational methods in a hybrid approach, which has found applications in quantum chemistry [131,132] and also lattice gauge theories [133].Particularly in the context of lattice field theories, where one eventually aims at recovering the continuum limit, increasing the size of these quantum simulators will be important in the near future.In addition, including both the matter particles and interaction carriers in the quantum simulation will also be important, allowing one to get closer to the higher-dimensional non-Abelian gauge theories of the standard model, which is one of the longer-term goals [134][135][136][137][138][139][140][141].
In the analog approach to quantum simulation [41][42][43][44][44][45][46][47][48][49][50][51], these spins are coupled to the transverse phonons of the ion crystal by applying an offresonant spin-dependent dipole force.This force can be obtained by a bichromatic laser beam with a pair of tones n = 1, 2 of frequency (wave-vector) ω L,n (k L,n ), which are either (i) symmetrically detuned with respect to the red and blue motional sidebands [38,39], or (ii) far-detuned from the atomic transition [142,143].In the following, we consider the second scheme although, as discussed later, the formalism also applies to the former, albeit in a different spin basis.
In contrast to the accumulation of errors in digital quantum simulators, which arise from both the imperfect operations in a gate sequence and the approximations inherent to the Suzuki-Trotter expansion, it is more difficult to account for the growth of errors in analog quantum simulators.Nonetheless, one expects that their accumulation in time will not be as detrimental as in a generic digital approach, especially when one is interested in recovering intensive observables [16].Accordingly, the common expectation is that one will be able to demonstrate quantum advantage using near-term experiments with these analog quantum simulators [20,21].In fact, some experiments have already been able to track real-time dynamics of a many-body model, going beyond the capabilities of current classical computers with stateof-the-art numerical algorithms [144].Even if it is difficult to provide mathematical proofs of quantum advantage, as one is departing from the quantumcomputing framework in which the scaling of required resources for a target accuracy is routinely estimated, there has been recent progress in this direction [145,146].
Working in the Lamb-Dicke regime [143], the lightmatter interaction of the ions with the bichromatic laser beam leads to a local interaction between the spins and the ion displacements.If the laser illuminates the entire ion chain, this reads where we have introduced the beatnote frequency (wave-vector) ), and the Pauli operator In the above expression, the force strength reads g = ℏΩ L ∆k L,z , where Ω L is the differential ac-Stark shift between the two electronic states [143].For simplicity, we will assume that ∆k L ||e z , such that ∆k L,x = 0 from now on.Working in the weak-force regime it is possible to obtain an effective spin model with long-range interactions mediated by the transverse phonons [40], which governs the slower dynamics of the spins [41][42][43][44][44][45][46][47][48][49][50][51].Using the coarse-graining in Eq. ( 16), this coupling can be expressed as which can be understood as a Yukawa-type coupling if one writes σ z i = 2ψ † (t, x i )ψ(t, x i ) − 1 in terms of a local fermionic field.Here, we have introduced the harmonic source terms Let us now address an important point by discussing when the coarse-grained description is expected to capture the properties of the long-range spin-spin interactions.The idea is that, whenever the harmonic sources (24) oscillate at a frequency that is close to the frequency ω(k) (8) of the lowest zigzag mode at k = π/d, namely ∆ω L ≈ ω(π/d), then the long-wavelength approximation will provide reliable results.This is actually more general than working at the vicinity of the structural phase transition, which is a low-energy approximation since the zigzag mode becomes the soft mode of the transition ω(π/d) ≈ 0. Accordingly, the long-wavelength approximation is also valid at other parameter regimes far from the structural transition, in which the additional non-linearities are unimportant.It is in these regimes, in which the elastic terms (3) suffice to describe the problem, where most of the experimental trapped-ion quantum simulators work [41][42][43][44][44][45][46][47][48][49][50][51].
In this case, the coarse-grained theory corresponds to a Klein-Gordon field with the dispersion relation of Eq. ( 8), and the real-time dynamics of the spins is governed by a unitary evolution operator with an effective Ising Hamiltonian where we have introduced the spin-spin coupling strengths As first noted in the experimental works [43,44], by controlling the detuning of the laser beams used to generate the long-range interactions in a trapped-ion crystal, the decay of the spin-spin couplings with the inter-ion distance can be approximately fitted to a power law with a tunable exponent.In the expression above, the distance decay of the spin-spin couplings is controlled by a dimensionally-reduced Euclidean propagator of the Klein-Gordon field, which is obtained after integrating the temporal components in x 1 − x 2 = (τ, x i − x j ), and reads as follows This propagator can be interpreted in terms of excitations of the scalar field with an effective mass/Compton wavelength where we recall that the trapped-ion bare mass has been defined in Eq. (13).The interpretation in terms of the coarse-grained picture is really transparent, the spin-spin interactions are mediated by the Klein-Gordon field, and the spin-spin couplings are controlled by the distance decay of the corresponding propagator and, thus, by the effective Compton wavelength.It is important to emphasize that, as will be discussed further below, a classical numerical simulation of the real-time dynamics (25) of the full Yukawacoupled QFT of spins and scalar fields is a very complicated problem especially in the presence of further non-linearities, finite temperatures, and certain microscopic details of the trapped-ion system.This connects to the prospects of using experimental quantum simulations in the context of quantum advantage, as mentioned in the introduction.
In the continuum QFT, the propagator (27) can be expressed in terms of modified Bessel functions [147] which, for the case of one spatial dimension, lead to exponentially-decaying spin-spin couplings with a typical decay length controlled by the Compton wavelength [53].For a standard lattice regularization of the Klein-Gordon QFT, which only has nearestneighbor couplings leading to Eq. ( 10), the previous integral can be evaluated by extending the momentum to the complex plane kd → z ∈ C, and by noticing that the integrand contains a simple pole that contributes with an exponential distance decay [147].However, for the trapped-ion regularization, the dispersion relation (8) also presents a branch cut, which contributes with an additional term with a dipolar distance decay.Altogether, the spin-spin couplings read where we have introduced an effective strength and η x = k L,x ℏ/2m a ω x is the Lamb-Dicke parameter.In the final expression (29), one can substitute the inhomogeneous equilibrium positions x i that stem from the solution of Eq. (2).Given the relation of the bare mass with the zigzag mode frequency in Eq. ( 13), we see that the effective Compton wavelength (28) that sets the range of the interactions in turn depends on by how close the laser beatnote is with respect to this vibrational mode.
As discussed in detail in [53], the exponential part of the spin-spin interactions is the typical contribution of a Yukawa interaction in D = 1 + 1 dimensions [147,148].On the other hand, the physical lattice regularization of this QFT stemming from the trapped-ion microscopic model also includes longrange couplings that modify the dispersion relation.This leads to a branch-cut discontinuity that is responsible for an additional dipolar-decaying part in the spin-spin couplings.As discussed in [53], the effective coarse-grained description of these Yukawa-type interactions is very accurate already for moderate-size chains with tens of ions, considering realistic parameters and inhomogeneous ion crystals.
Since the distance dependence of the spin-spin couplings can be inferred from various experimental techniques that measure the real-time dynamics of the effective spin model [43,45,46], it follows that one could use the spins as probes of the underlying effective relativistic QFT.One can thus use experiments in which the crystal sizes are sufficiently large to admit a continuum limit as quantum simulators of a relativistic Yukawa-type problem.In analogy to the simulations of lattice gauge theories [33][34][35] mentioned above, the trapped-ion spins would mimic the matter degrees of freedom since, under a Jordan-Wigner transformation [149], they can be interpreted as fermionic matter.In the present case, this matter would be quenched, and sit on the sites of a real physical lattice corresponding to the ion crystal.Instead of having a gauge field to mediate the interactions, which must be defined on the links of the lattice to allow for a local symmetry, here one has a global inversion symmetry of a real scalar field ϕ(x) defined on the lattice sites.Note that the continuum limit is not recovered by sending the physical lattice spacing to zero, but rather by working in parameter regimes where only the long-wavelength properties are of relevance.These regimes correspond to the vicinity of a second-order phase transition.
To test the validity of this expression (29) in a specific realistic setup, one may consider crystals of up to N = 51 atomic 40 Ca + ions in a linear Paul trap, forming a chain with an overall length of L ≈ 250 µm [150].Reference [150] reviews several aspects of the coherent manipulation of such long strings.It is also shown that stability of the chain is characterised by a lifetime of approximately 27 s, a period after which the collisions of the ions with the background gas become important and can even melt the crystal [151,152].This sets an upper bound for the possible time of the experimental runs, and demands fast and efficient laser cooling in order to reach a low vibrational state of all axial and transverse modes, as required for highfidelity quantum control.We note that these constraints can be met by using resolved-sideband cooling for the transverse modes [153], and polarizationgradient cooling for the axial ones [154].For the transverse modes relevant for our work, we will assume that the mean phonon number ranges between 1-30, which will become important below.
The spins are encoded as optical qubits in the long-lived electronic states ↓ = 4S 1/2 (m = +1/2), and ↑ = 4D 5/2 (m = +5/2), which have long coherence T 2 = 64 ms and decay T 1 = 1 s times [150].The motional degrees of freedom depend on the average lattice spacing d ≈ 5 µm, and the transverse and axial trap frequencies, which have typical values of ω y /2π = 2.93 MHz, ω z /2π = 2.89 MHz and ω x /2π = 127 kHz in most of the experiments discussed in [150].In the following, we assume a N 1 = 30 ion chain, in which the classical estimate of the linearto-zigzag transition (15) for a fixed axial trap frequency of ω x /2π = 127 kHz corresponds to the critical transverse frequency ω z,c /2π ≈ 2.7 MHz.We thus consider modifying this trap frequency in the range ω z /2π ∈ [2.75, 2.89] MHz to cover the regimes in the symmetry-preserved phase, i.e. a stable chain configuration, which are well described by either the Klein-Gordon effective QFT for ω z ≫ ω z,c , or otherwise by the full λϕ 4 model for ω z ≳ ω z,c .
In the experiments discussed in [150], they use a bichromatic laser scheme with a pair of beams symmetrically detuned with respect to the red and blue motional sidebands [38,39].This leads to a state-dependent dipole force that acts in a different spin-basis with respect to Eq. (21).As discussed in [40,150,155], the spin-spin interactions mediated by the transverse phonons are given by the following expression where the recoil energy is E R = (ℏ∆k L ) 2 /2m a , and the laser beatnote is now referenced to the qubit transition ω L,1 − ω L,2 = ω 0 + ∆ω L , with ω 0 /2π = 411.5 THz.Also, in contrast to the previous case (21), Ω L is now the quadrupole Rabi frequency [156] instead of the differential ac-Stark shift.In this expression (31), we have introduced the normal-mode frequencies ω z,n1 and wavevectors M in1 of the transverse phonons [70].We note that Eq. (31) does not rely on the approximations used to obtain the coarsegrained QFT prediction (29).The only common implicit assumption is that one neglects the higher-order quartic terms, as well as other off-resonant contributions beyond the dipole force (21).
Considering the experimental trap frequencies ω z /2π = 2.89 MHz and ω x /2π = 127 kHz, and setting ω L,1 − ω L,2 = ω 0 + ∆ω L where ∆ω L /2π = 2.43 MHz is red-detuned with respect to the zigzag mode at ω z,N1 /2π = 2.7 MHz, we compare in Fig. 2 the spin interactions in a trapped-ion chain of N 1 = 30 ions obtained from the exact expressions of the inhomogeneous crystal (31) with those of the coarse-grained Yukawa-mediated expression (29).In this figure, we have set Ω L /2π = 1 MHz, such that J i0,i0+1 /h ≈ 1.4 kHz for the ion at the center of the chain i 0 = 15.The agreement of the coarse-grained QFT prediction, which has no fitting parameter, is rather remarkable for such moderate-size chains.As expected, there are some deviations at distances around the UV lattice cutoff, but the accuracy becomes very good as the ion distance increases.
In summary, we can conclude that trapped-ion quantum simulators of spin models, like those of the experiments in [41,42,44,[44][45][46][47][48][49][50][51] but performed on chains with a few tens of ions and a spin-dependent force that is red-detuned with respect to the transverse zigzag mode, can be used to probe the physics of a relativistic QFT when interpreted in the light of a Yukawa-type interaction.As seen from the glasses of QFT, the problem becomes more interesting in the presence of quartic interactions and non-zero temperatures, as discussed below.

Effect of λϕ 4 term and non-zero temperatures
Let us now discuss how to take into account the nonlinearities (6) that go beyond the harmonic approximation (3), as well as a non-zero temperature, both of which become relevant in a trapped-ion experiment.In the parameter regime underlying Fig. 2, the trapping conditions are far from the linear-to-zigzag critical point ω z ≫ ω z,c , such that the harmonic-crystal approximation (3) is an accurate description of the collective phonons.In terms of the coarse-grained QFT (6), the bare quartic term (20) is negligible in comparison to the bare mass (13), and the effective QFT can be reduced to that of a real Klein-Gordon field.Temperature enters in the conditions that must be imposed on the strength of the Yukawa-type coupling (23), which we recall had to fulfil the weakcoupling constraint (22) in the zero-temperature limit.In the presence of thermal fluctuations, there can be some bosonic enhancement of the Yukawa-type coupling, and the condition (22) must be upgraded to where n(k) is the mean excitation number of the scalar field.In this regime, the real-time timeevolution of the spins is still described by a longrange Ising model (25) in complete analogy to the zero-temperature case.Once again, the spin-spin couplings (26) are proportional to a dimensionallyreduced Euclidean propagator (27).The important difference is that this is not the free propagator of a Klein-Gordon QFT, and can get a temperaturedependent contribution as one goes beyond the harmonic non-interacting limit.
As mentioned previously, one can control n(k) by means of laser cooling.For a single trapped ion [157], one can show that the resulting state corresponds to a thermal Gibbs state, and that the mean-phonon number can be controlled by adjusting the ratio of the cooling and heating rates, which depend themselves on the detuning and intensity of the cooling laser.For a trapped-ion chain, the steady state that results from the cooling will depend on the laser-cooling scheme.For instance, in the resolved-sideband limit, one can individually cool each of the transverse modes to a target mean excitation number n(k).In this work, we assume that, after such cooling stage, the trappedion crystal is allowed to thermalize, and one can then define a single effective temperature according to the Bose-Einstein distribution of a thermal Gibbs state Given the nature of the coarse-grained approximation underlying (6), we believe that, even if there are deviations to this idealised thermalization, and the effective temperature varies for the different modes, the only relevant quantity is the effective temperature around the zigzag mode k = π/d.Let us now discuss why we need to go beyond the harmonic limit.In this limit, increasing the temperature only results in more stringent conditions for the Yukawa-coupling strength (32), which in turn result in a weaker magnitude for the spin-spin interactions (29).The situation becomes more interesting in the presence of non-linearities, such as the quartic coupling (20), which become more important as one approaches the linear-to-zigzag transition ω z → ω z,c .At zero temperature [53], a path integral can be used to show that the spin-spin couplings (26) are still described by Eq. ( 26), but this time controlled by a dressed propagator and renormalised sources.This propagator depends on all of the possible scattering processes of the self-interacting scalar field, as its excitations propagate between a pair of distant ions.These interaction effects will shift the physical mass of the carrier from m 2 0 → m 2 P , changing the effective Compton wavelength (67), which will in turn change the range of the Yukawa-type interactions (29).This can be inferred experimentally from changes of Fig. 2 as one approaches the structural transition.
Perturbatively, one expects a zero-temperature shift of the bare mass that scales with the quartic coupling and stems from the so-called tadpole Feynman diagram δm 2 0 | T =0 ∝ λ 0 [147,148], which will be discussed at length below.As estimated in [53], these type of effects are inhibited by the very large rigidity of the ion chain, as the quartic coupling scales with the inverse fourth power of the rigidity modulus (20).Moreover, given the fact that the effective coarsegrained parameters in Eqs.(11), ( 13), (17) and (20) depend on the microscopic experimental parameters in a convoluted manner, it is not straightforward to modify λ 0 independently of the others to see its effect on the range of the interactions.
In addition to the above zero-temperature shift of the bare mass, the tadpole diagram also contributes to the so-called thermal mass [158,159].Perturbatively, this reads as follows where the proportionality stems from changes that would be required to convert from natural units into SI units.Regardless of the proportionality factor, the important result of thermal field theory is that, in spite of having a small quartic coupling, the above integral can lead to a shift proportional to some power of the ratio δm 2 0 | T ∝ λ 0 (T /m 0 ) α [158].For D = 1 + 1 dimensions, we find δm 2 0 | T ∝ λ 0 (T /|m 0 |) in the hightemperature regime T ≫ m 0 , which can thus amount to a large shift even for a perturbative λ 0 .In the present context, this thermal shift will change the Compton wavelength (28), and the range of the spinspin interactions (29); an effect that will be characterised non-perturbatively below.
Coming back to the trapped-ion regularization of the λϕ 4 QFT (6), the classical critical point (15) of the Z 2 -breaking phase transition, obtained by setting the bare mass (13) to zero m 2 0 = 0, will get contributions from the thermal masses.Hence, the physical mass will be dressed with both temperature and quartic-coupling contributions, such that the phase transition at m 2 P (λ 0 , T )| c = 0 now corresponds to a critical surface in parameter space (m 2 0 , λ 0 , T ).Determining how the critical point flows with temperature and coupling strength in the lattice model is a non-perturbative problem that requires going beyond the previous discussion and will be addressed in the following sections.In general, if one starts in a symmetry-broken phase m 2 P (λ 0 , T ) < 0, by solely increasing the temperature, symmetry restoration can take place m 2 P (λ 0 , T + δT ) > 0 [63,64], such that one ends in a symmetry-preserved phase.In terms of a linear-to-zigzag phase transition at finite temperatures [160], an analogue of this restoration of symmetry has actually been observed already in experiments with small trapped-ion chains [161][162][163].To the best of our knowledge, the connection to relavistic QFTs and thermal masses has not been previously noticed in the trapped-ion literature.In these experiments, the cooling lasers are not only used to prepare an initial thermal state, but are actually applied continuously during the experiment, such that one is exploring the steady state of a driven-dissipative system.In spite of these differences, the observations resemble the phenomenon of thermal masses and the restoration of symmetry.As discussed in [161], a linear shift of the critical ratio κ z,c (15) with temperature has been reported, which is somewhat reminiscent of the previous thermal mass shift.In the following section, we will present a self-consistent non-perturbative method that can be used to derive quantitative predictions of how the critical point flows, and how the range of the spin-spin interactions changes with temperature.

Self-consistency beyond mean field theory
In this section, we present a detailed account of our self-consistent prediction for the range of the spinspin interactions, and how it changes with temperature as one approaches the linear-to-zigzag transition.We start by reviewing the functional approach to the self-interacting scalar field theory, and then move on to discuss our approach to get a set of finite self-consistent equations in spite of UV and IR divergences.

Perturbative generating functional of λϕ 4 fields
In this subsection, we review the functional approach to the diagrammatic perturbation theory of the selfinteracting scalar field [147,148,164].The central object in this approach is the generating functional, obtained by adding a source term to the path integral.In its Euclidean version, where the x = (τ, x) is obtained from x = (t, x) after a Wick rotation τ = it, the generating functional is given by where L the Lagrangian density associated to the Hamiltonian field theory in Eq. ( 5), provided one works in imaginary time.The generating functional gives the n-point Green's functions upon derivation with respect to the sources where the expectation value is taken on the vacuum, and is the normalized generating functional.
In the absence of interactions, the Lagrangian L = L 0 reduces to a real Klein-Gordon theory quadratic in the fields, and the normalised generating functional finds a simple analytical expression where ∆ 0 (x − y) = ⟨ϕ(x)ϕ(y)⟩ is the Euclidean propagator of a free scalar field.The propagator is most conveniently written in terms of its Fourier decomposition where the Euclidean momentum k = (k 0 , k) is related to the 2-momentum k = (ω, k) in Minkowski spacetime by k 0 = −iω.In an interacting theory L = L 0 + L int , such as L int = λ0 4! ϕ 4 in our case, the generating functional can be obtained using the following identity in functional analysis This expression can be expanded in a power series in λ 0 to the desired perturbative order.The expansion can be graphically represented in terms of Feynman diagrams.Taking into account that the denominator in Eq. ( 39) cancels out the so-called vacuum diagrams, i.e. those diagrams that do not contain any external source terms, the final expression reads The crosses × = J(x) represent sources, the blobs = λ 0 interaction vertices, and the different free propagators obey = ∆ 0 (x − z 1 ), = ∆ 0 (0) and = ∆ 0 (z 1 − z 2 ).As usual, we recall that one has to integrate over the location of the interaction vertices, here denoted by {z i }.

Feynman diagrams and self-consistent equations
The Feynman diagrams in Eq. ( 40) provide quantum corrections to the n-point correlation functions of the QFT.For instance, the propagator of the interacting theory, which corresponds to the 2-point function, can be written as where Σ(k), the self-energy [165][166][167], contains the contributions of all 1-particle irreducible diagrams with two external legs.We recall that these diagrams are those that cannot be separated into disconnected pieces by cutting an internal propagator [148].At order λ 2 0 , the self-energy is given by the first, second and fourth diagrams in Eq. (40).The third diagram is not 1-particle irreducible and, thus, does not contribute to Σ(k).The first two diagrams belong to the tadpole type, namely, each loop integral depends on a single internal momentum that is independent of the external momentum of the propagator.As a consequence, the tadpole amounts to a renormalization of the mass and, in the (1 + 1)-dimensional case, is responsible for the only ultra-violet divergence of the QFT.The fourth diagram is the so-called sunrise or sunset diagram and, in contrast, also depends also on external momenta.In this way, this diagram contributes to both a mass and a wavefunction renormalization.
Evaluating the self-energy non-perturbatively is of course out of reach.It is however possible to resum the tadpole family to all orders in the quartic coupling.The result is encoded in the self-consistent equation This approximation to the self-energy becomes exact for an O(N ) vector model in the limit N → ∞ [168].
In condensed matter, it corresponds to the Hartree method of mean-field theory [169].At this point, we must discuss the occurrence of divergences in the QFT (6).In fact, the integral in Eq. (42) contains a UV logarithmic divergence in D = 1 + 1 dimensions, which must be regularized by introducing a cutoff.Since we ultimately want to use this selfconsistent resumation of Feynman diagrams to predict changes of the Yukawa-type spin-spin interactions in the trapped-ion chain (29), the regularisation scheme should correspond to a lattice.In the standard approach of lattice field theories, where the spatial derivatives are exchanged for discrete differences (9) that only lead to nearest-neighbor couplings, the continuum propagator appearing in Eq. (41) with the tadpole-resummed self-energy (42), must be substituted by where the analogue of the spatial momentum is As already mentiond above, the quasi-momentum lies within the First Brillouin zone k = 2π N1a n 1 for n 1 ∈ {1, • • • , N 1 }.In the propagator, we have also defined the tadpole-renormalised mass through the following self-consistent equation, sometimes refereed to as the gap equation, We note that the integrals over quasi-momenta are to be understood as mode sums dk 2π → 1 N1a n1 .In the thermodynamic limit, one sends N 1 → ∞, such that k ∈ [0, 2π a ), and the corresponding integral can be evaluated analytically Here, K(x) = π/2 0 dθ(1 − x sin 2 θ) −1/2 is the complete elliptic integral of the first kind [170], which is finite, showing that the UV divergence is regularised by the non-zero lattice spacing.For our scalar field theory, this is the only UV divergence.In Sec. 4, we will discuss how this self-consistent equation, as well as the expression that follow, can be adapted to the trappedion case where the dispersion relation (8) includes the not only nearest-neighbors but the full dipolar tail.For the moment, however, we continue with the standard Hamiltonian lattice regularization, and the corresponding propagator (43).
Let us now discuss the connection of this renormalised mass with the aforementioned phase transition.The physical mass of the interacting QFT is determined by the pole of the quantum corrected propagator (41) in Minkowski spacetime.The analytical prolongation is simply achieved by replacing k 0 → −iω, such that the propagator in Minkowski spacetime ∆(k) is obtained from −i ∆(−iω, k) → ∆(k).In the lattice version (43), the Lorenz invariance between energy and spatial momentum is broken.The physical mass is then defined as the on-shell energy at vanishing spatial momentum m 2 P = ω2 , and determined by the pole of the propagator In the mean-field tadpole approximation, the physical mass would thus be m 2 P = µ 2 , which implies that the classical critical point for the Z 2 -breaking phase transition m 2 0 = 0 will flow to a different value that can be obtained by solving for µ 2 = 0.When trying to solve this equation, one faces a different type of divergence in the QFT, namely an infra-red (IR) divergence.Note that the elliptic integral in Eq. (46) inherits the logarithmic IR divergence of the tadpole (42) when µ 2 = 0, since lim x→1 K(x) = ∞.This prevents criticality to be achieved for any finite negative value of the bare mass m 2 0 , independently of the value of the quartic coupling λ 0 and the non-zero lattice spacing.Accordingly, the vacuum of the 1+1 λϕ 4 theory would always remain in the unbroken phase if one sticks to this mean-field tadpole approximation, which is clearly wrong in light of other studies [88][89][90][91][92][93][94][95][96].We note that this caveat only appears in (1 + 1) dimensions, and forces us to go beyond the mean-field approximation.
In order to circumvent this problem, and obtain a line of critical points in parameter space (m 2 0 , λ 0 ), further quantum corrections need to be included in the self-energy.We will make the self-energy exact at second order in the quartic coupling by adding the sunrise contribution, which has a diagrammatic representation given by the fourth Feynman diagram in Eq. (40).Contrary to the previous tadpole terms, this contribution depends on the external momenta of the propagator, making a full self-consistent treatment that includes tadpole-and sunrise-like diagrams to all orders of the quartic coupling impractical.We can, however, include all tadpole decorations in the internal propagators of the sunrise diagram, leading in this way to an improved self-energy Σ(−iω, k).This can be achieved by introducing again the tadpoleresummed propagator in Eq. ( 43) in the mentioned self-energy, including thus the tadpole-like decorations to every order of the quartic coupling.Paralleling our discussion around Eq. (47), the critical line determined by m 2 P = 0 will thus be defined by the new condition where tadpole corrections come from Eq. (45), and d 2 kd 2 q (2π) 4 ∆td (k) ∆td (q) ∆td (k + q).(49) Eventually, we will also be interested in moving out of criticality, since it is the non-zero value of the physical mass m 2 P via its associated effective Compton wavelength ξ eff,P , which controls the range of the spin-spin interactions (29).Considering that the effective λϕ 4 model is only relevant close to the structural phase transition of the ion chain, we can assume that in the region of interest for the experiment, the physical mass m 2 P will be small.In order to solve for the pole of the improved propagator in Eq. (47), it will then suffice to consider where the Taylor expansion of the sunrise diagram yields Close to this pole, the propagator with contributions from all tadpoles and the sunrise diagram has the expression where we have rotated back to Minkowski spacetime.
Here, one can readily identify the physical mass to be In addition to the additive contribution of these Feynman diagrams to the physical mass, one also observes the appearance of a multiplicative contribution from the residue at the pole We thus see that quantum corrections do not only modify the mass of the theory, but also the normalization of the field itself.This is the physical meaning of z, the so-called wave-function renormalization.In order to have the field canonically normalized, it is necessary to rescale ϕ(x) → √ zϕ(x).This explains the appearance of z in the physical mass (53).
Let us now discuss the improved prediction of the critical line, which is obtained by solving m 2 P = 0.This is no longer impeded by the IR divergence of the tadpole integrals, as the resummed contributions to the mass µ 2 are no longer required to be zero.In some sense, the lattice regularization provides a cutoff at large momenta that allows us to get finite results in the UV limit, while the addition of the sunrise diagram provides an effective "mass cutoff" at low energies that allows us to get finite results in the IR limit.In summary, the equation that must be solved is given by where the wavefunction renormalization z is given by Eqs. ( 54) and (51), and the tadpole contribution µ 2 is that of Eqs. ( 45) and (42).Even though Eq. ( 55) is not a proper self-consistency equation when written in this form, we opt for simplicity and refer to the whole set of equations ( 45) and (55) as the self-consistency equations.In the following sections, we will present numerical solutions of this set of equations applied to the Yukawa-type interactions in trapped ions.Let us, however, first discuss how the formalism of thermal field theories can account for non-zero temperatures.

Non-zero temperature and thermal field theories
As stated in Subsec.2.3, the main goal of this work is to explore the effect of non-zero temperatures in the Yukawa-mediated spin-spin interactions of a trappedion quantum simulator (25).We argued that, in the vicinity of a structural phase transition, the nonlinearities will lead to a shift of the physical mass of the effective λϕ 4 QFT, which will change the distance decay of the spin-spin couplings (29).Moreover, there can be additional contributions at non-zero temperatures related to the perturbative thermal mass of Eq. (34).In this subsection, we discuss how to generalise the previous self-consistency equations (55) to a non-zero temperature using the formalism of thermal field theories [158,159].
According to a path integral approach [147] in Euclidean time for a system at thermal equilibrium [171], the generating functional of Eq. ( 35) corresponds to the partition function of the λϕ 4 model in the limit of a vanishing temperature.For non-zero temperatures T > 0, the integrals in Eq. ( 37) must be modified, as the spacetime x = (τ, x) is no longer the Euclidean plane.The field in the path integral must fulfill periodic boundary conditions where β = 1/T is the inverse temperature in natural units, such that the Euclidean plane is effectively compactified into a cylinder R 2 → S r × R of radius r = β/2π.Given the periodicity of the scalar field (56), the propagator obeys the Kubo-Martin-Schwinger condition [172,173], namely ∆ 0 (τ, x) = ∆ 0 (τ + β, x).Accordingly, the transformation to momentum space (38) requires using a Fourier series [158,159,174] in the temporal coordinate instead of a Fourier transform where ω n0 = 2πn 0 /β are the bosonic Matsubara frequencies.The important aspect of the Matsubara formalism is that the generating functional of the equilibrium n-point functions of the thermal field theory has the same functional form as the T = 0 case (39), provided one substitutes the frequency integrals by a series in the Matsubara frequencies.Accordingly, one can use the previous diagrammatic results (40), as well as the self-consistency equations (55) with the lattice propagator (43), by substituting: , where n 0 ∈ Z and n 1 ∈ {1, • • • , N 1 }; (ii) the Euclidean lattice propagator by the Matsubara propagator ∆ 0 (k) → ∆ 0 (ω n , k); and (iii) the momentum integrals by mode sums Let us illustrate this procedure by considering the lattice-regularised tadpole contribution (42).At finite temperatures, we have  0) is proportional to a dimensionless ratio µ 2 Σsr(0)/λ 2 0 that is, again, independent of the quartic coupling.We compare the exact T = 0 result in the thermodynamic limit (black dotted line), discussed in Appendix B, with the finite-temperature expression obtained by performing the Matsubara sums explicitly (73) (see Appendix A), for a lattice with N1 = 30 sites (red dashed line).
The Matsubara sum can be performed explicitly, yielding This contribution is the mean field counterpart of Eq. ( 34), as detailed in Appendix A. We plot in Fig. 3(a) the dimensionless quotient Σ td /λ 0 for N 1 = 30 as a function of the temperature, measured in lattice units, and for the value of the tadpole renormalized mass µa = 1.The dotted line in that figure corresponds to the zero temperature thermodynamic limit, evaluated in Eq. (46).We observe a quite small deviation from the thermodynamic limit already for the moderate number of sites we have chosen.
Remarkably, the Matsubara mode sums involved in the sunrise contribution Σ sr (0) (53) can also be performed analytically.The resulting expression can be found in Eq. ( 73) of Appendix A. The dimensionless combination µ 2 Σ sr (0)/λ 2 0 is shown in Fig. 3(b), where again a small deviation from the zero temperature thermodynamic limit is found for N 1 = 30.Conversely, the wavefunction renormalization (54) at finite temperature needs to be computed numerically, as no analytical expression for the Matsubara mode sums was found.This requires the truncation of the Matsubara sums to a finite number of modes n 0 ∈ {−N 0 , • • • , N 0 }, and thus raises the issue of the convergence with N 0 .This is studied in Appendix A, obtaining a fast convergence which will be important for the efficiency of our numerical approach, even in the typically more demanding low temperature limit.

Critical line and trapped-ion spin-spin couplings
In this section, we will numerically solve the selfconsistent equations (45) and (55) at non-zero temperatures.We will obtain results for the standard Hamiltonian lattice discretization of the scalar field in Eqs (43)-( 44), but also introduce the dipolar tail of the dispersion (8) to be able to make explicit predictions for the trapped-ion case.

Numerical estimate of the critical line
As discussed above, the Z 2 -breaking phase transition in the λϕ 4 model is characterised by a classical critical point at m 2 0 = 0, which will flow with temperature and quartic coupling due to thermal and quantum effects that shift the pole of the propagator to the physical mass m 2 0 → m 2 P .As a consequence, the classical critical point will become a critical surface in parameter space (m 2 0 a 2 , λ 0 a 2 , T a) determined by solving the equation m 2 P (λ 0 a 2 , T a) = 0.If one fixes the temperature to a specific value, we want to obtain the corresponding critical line for the bare mass m 2 0 | c as a function of the bare coupling strength λ 0 .This line will separate the symmetry-broken m 2 0 < m 2 0 | c from the symmetry-preserved m 2 0 > m 2 0 | c phase at the given temperature T a.
In order to numerically obtain these critical lines for various non-zero temperatures, one needs to impose the condition m 2 P = 0 in Eq. ( 55), and solve the selfconsistency equations to extract the bare critical parameters (m 2 0 a 2 , λ 0 a 2 )| c .We note that the wavefunction renormalization, contributing multiplicatively to the physical mass (53), does not play any role in the determination of the critical points.Then, the routine followed to compute the critical line reads as follows Routine R1 1 Select an interval µ 2 1 a 2 , µ 2 2 a 2 in the µ 2 a 2 axis 2 For each µ 2 a 2 in µ 2 1 a 2 , µ 2 2 a 2 and a given T a: 3 Impose m 2 P = 0 in the lattice counterpart of Eq. ( 55) using Matsubara sums to obtain λ 0 a 2 as a function of µa, T a and N 1 4 Substitute µ 2 a 2 in Eq. ( 45) to obtain In Fig. 4, we represent the critical line for different temperatures.For a fixed value of the temperature T a, the region above the critical line represents a symmetric phase, adiabatically connected to the Klein-Gordon thermal state, whereas the region below is the symmetry-broken phase in which the scalar field acquires a non-zero expectation value and can no longer be adiabatically connected to the Klein-Gordon vacuum.As can be observed in this figure, the critical lines take lower values of the quartic coupling as the temperature increases.This behaviour is consistent with the appearance of a thermal mass (34) and the phenomenon of restoration of symmetry [63,64].For a fixed value of the coupling strength λ 0 a 2 , if one starts at a point (m 2 0 a 2 , T 1 a) in which the symmetry is broken (i.e. a point below the corresponding critical line of Fig. 4), one may increase the temperature towards T 2 > T 1 , such that the corresponding equilibrium state lies now above the critical line, and thus belongs to the symmetry-preserved phase.Let us note that the dependence of the critical line with the quartic coupling is similar to the zero-temperature results obtained using other numerical methods, such as Monte Carlo [88,91].If one aims at exploring regimes beyond those of relevance in the trapped-ion realization (see our discussion in Sec.2.3), there are certain limitations of the current approach that are discussed in detail in Appendix C.
Once the critical lines have been derived, a more  55) for the critical point m 2 P |c = 0. We consider a lattice regularization (43) with N1 = 30 sites and use the analytical Matsubara mode sums discussed in Appendix A along with the Routine R1 described in the text.When the temperature increases, the effect of a thermal mass makes the critical line take lower values for the quartic coupling.For each colored line, the upper region corresponds to the symmetry-preserved ground state (SP-vacuum), whereas the region below corresponds to the ground state with spontaneous breaking of the parity symmetry (SSB-vacuum).
informative figure for the connection to the trappedion case and the Yukawa-mediated interactions discussed in the following section would be a contour plot for the non-zero values of the physical mass m 2 P a 2 in the plane (m 2 0 a 2 , T a) at a fixed quartic coupling λ 0 a 2 .This will be the first step to make a connection with the distance decay of the spin-spin interactions in future sections.Achieving this goal requires incorporating the wavefunction renormalization (54), which contributes multiplicatively to the physical mass (55), into a different numerical routine Routine R2 1 Select the interval of dimensionless temperatures (T 1 a, T 2 a) and a fixed λ 0 a 2 value 2 For each T a in (T 1 a, T 2 a): Update m 2 P a 2 z −1 and compute m 2 0 a 2 from Eq. (55) and Eq. ( 45) respectively 7 Compute z using (54) In practice, the same ϵ variation of µ leads to different size increments on m 2 P a 2 and m 2 0 a 2 , depending on T a and the previous m 2 P a 2 value.To obtain uniform increments independently of position in parameter space, we heuristically define an adapted step ϵ i ∝ m 2 P,i−1 a 2 /(1 + T a) (1/3) .Now, we can run a simulation based on this numerical routine to obtain a contour plot of the physical mass.The numerical results are shown in Fig. 5, where the region in parameter space with a non-zero physical mass is coloured according to the scale specified in the rightmost inset.The white region corresponds to the symmetry-broken phase, where the physical mass would become negative, leading to a transition from a single to a double well in the effective potential [85,86].This potential, via the self-consistency equations, has quantum corrections stemming from all the tadpole diagrams, as well as perturbative corrections of the sunrise diagram.We can see how the symmetry-broken region changes with temperature for the various quartic couplings.In this contour plots, the restoration of symmetry by increasing temperature becomes very clear, as it corresponds to any vertical line connecting the symmetry-broken and symmetry-preserved phases.
As explained before, our self-consistent equations are exact to order λ 2 0 , but miss many higher-order contributions.It is therefore important to benchmark the performance of our approach.We will consider a specific universal quantity, namely the ratio f c = λ 0 /µ 2 | c at the T = 0 critical point of the λϕ 4 QFT in the continuum, which has also been computed with other numerical approaches.In the continuum, the UV divergence of the tadpole contribution (42) would make the dimensionless ratio λ 0 /m 2 0 | c vanish.It is then customary to replace the bare mass by the UV-finite tadpole renormalised mass µ 2 (45).Using our lattice discretization, we can access the critical ratio as a function of the lattice spacing f c (a), or, alternatively, as a function of f c (λ 0 a 2 ).Its value in the continuum QFT can be thus obtained by estimating f c (λ 0 a 2 ) as the dimensionless coupling strength λ 0 a 2 → 0.
To compute f c , one can make use of Routine R1 with some minor modifications.For a fixed number of spatial lattice sites N 1 , and using a small-enough interval of values for µ 2 a 2 , one can naively extract the critical ratio by plotting λ 0 a 2 | c vs µ 2 a 2 | c , and fitting the graph to a linear function in order to extract the slope.We find f c ≃ 20.105 using N 1 = 2000.As discussed in detail in [175], logarithmic contributions must be taken into account in order to be get a more precise estimate.Guided by this work, we employ three logarithmic models to fit our results, as depicted in Fig. 6.As shown in this figure, the three models give the same critical ratio value up to the second decimal place, f c = 20.10.We again use N 1 = 2000, as we appreciate a clear convergence to the continuum (see the inset in Fig. 6).Within our approach, the universal ratio can actually be computed directly in the continuum.From equations (48) and (49), we have f c = 6(2π) 4 /I with In Appendix B we evaluate analytically this integral with the help of Feynman parameters obtaining f c = 20.1055,which certifies the good behaviour of the continuum limit of the lattice discretization.
The quantum critical point of the λϕ 4 theory in 1+1 dimensions has been studied with a variety of nonperturbative methods, such as Tensor Networks [175][176][177] or Monte Carlo [94], among others [178,179].All these works find values near f c ≃ 66.0 with increasing levels of precision.The deviation of our self-consistent approach from these more-accurate estimates is common to in quantum many-body models with critical points, as illustrated for instance by the mean-field underestimation J/h| c = 1/2 of the critical point h/J| c = 1 of the Ising model in a transverse field [180].Common to this type of mean-field treatments, our treatment overestimates the stability of the ordered symmetry-broken phase under an increase in the mass.We can say that the simplicity of our method, which reduces drastically the computational complexity to a solution of self-consistency equations, comes at a price.On the other hand, the interest of our approach lies in being the most economic one capable of coping with the divergences and still detecting the phase transition and study the phase di- agram.Indeed, this simplicity will be crucial in the next section, where we will apply it to the trapped-ion quantum simulators of spin models [41,42,44,[44][45][46][47][48][49][50][51].Indeed, the non-standard dispersion relation (8), responsible for the dipolar tail of the spin-spin interactions in Eq. ( 29), might be difficult to treat accurately with more precise methods, such as the aforementioned tensor networks which can approximate long-range couplings by using the so-called matrixproduct operators.Moreover, since we are ultimately interested in dealing with non-zero temperatures, this would further increase the computational complexity of tensor-network methods, requiring again to work with matrix-product operators.Finally, in the context of the trapped-ion implementation, addressing the complete problem of the Yukawa-type spinspin interactions would require solving the real-time dynamics of the spins under locally coupled to the phonons, which evolve on a much slower timescale in comparison to the effective λϕ 4 model.This would require simulation long-time dynamics, which is a notably difficult for Tensor Networks, and out of reach for Monte Carlo methods due to the sign problem.Our self-consistent treatment should thus be seen as a proof of concept for the proposal to use thermal effects to study interacting QFTs with trapped-ion simulators.

Estimates for the trapped-ion quantum simulator
Let us now discuss how the previous results can be connected to the trapped-ion quantum simulators of spin models [41,42,44,[44][45][46][47][48][49][50][51].In Sec. 2, we de-scribed in detail how, far from the linear-to-zigzag structural phase transition (15), the spin-spin couplings mediated by the transverse phonons (31) of a trapped-ion chain can be described accurately by the Yukawa-type interactions mediated by a Klein-Gordon field (29) (see the comparison in Fig. 2).As the trap frequencies are modified and one gets closer to the structural phase transition, the Coulomb nonlinearities start to play a bigger role, and transform this QFT into the λϕ 4 model (6) with an effective speed of light in Eq. (11), and bare parameters in Eqs.(13), (17) and- (20).We argued that the specific dispersion relation ω(k) for the transverse modes of the trapped-ion chain (8) is very similar to that of a discretized scalar field (10), which underlies the Feynman propagator (43) we used in the solution of the self-consistent equations (55) of the previous subsection.Note that this propagator, as well as the results in Figs.3-6, have all been obtained using natural units ℏ = c = k B = 1.Moreover, in the lattice discretization at finite temperature, these self-consistency equations (45) and (55) are rewritten in terms of finite mode sums, and depend on dimensionless parameters obtained through a specific power of the lattice constant m 2 0 a 2 , λ 0 a 2 , T a, In the trapped-ion case, the effective speed of light c t (11) only appears after the gradient expansion leading to Eq. ( 6).This gradient expansion cannot account for the branch-cut discontinuity of the dispersion relation (8) when extended to the complex plane, and would thus miss the dipolar part (29) of the Yukawa-mediated interactions.Therefore, rather than setting ℏ = c t = k B = 1 in the coarse-grained trapped-ion case, it would be better to work with the full phonon propagator prior to the long-wavelength approximation.This requires reformulating the selfconsistency equations (45) and (55) in terms of dimensionless trapped-ion parameters, which can no longer be obtained by multiplying the microscopic parameters with a power of the ion lattice spacing d, as we need to use SI units.
In order to find such a formulation, we start by revisiting the tadpole-resummed propagator on the lattice (43).For trapped ions, the tadpole-resummed propagator is analogous to Eq. ( 43), but has inverse squared-energy dimension where k 0 has dimensions of inverse time, and will be substituted by the Matsubara frequencies for a nonzero temperature.Therefore, the trapped-ion analogue of the spatial lattice momentum k in Eq. (44) must also have units of inverse time.Moreover, in the absence of quartic interactions, the effective bare mass in Eq. ( 13) should appear as a pole in this propagator µ 2 = m 2 0 upon the substitution of k 0 → −iω.Taking into account these two conditions, we find that the analogue of the lattice spatial momentum (44) that appears for the nearest-neighbor discretization of the scalar field (9) (62) where we have used the truncated Riemann zeta in Eq. (14).Accordingly, k has the desired dimension of inverse time, and one sees how the dipolar tail of the phonons dispersion relation (8) enters in the propagator.
Once the propagator has been identified, we can rescale the tadpole self-consistency equation (45) with a certain power of the effective speed of light and Planck's constant, such that the equation has the right dimension of energy squared.In SI units, the Matsubara frequencies in Eq. (57) are expressed in terms of ω n0 = 2πn 0 /βℏ, where β = 1/k B T , such that the trapped-ion tadpole self-consistency equation reads (63) Considering that the quartic coupling (20) in SI units has dimensions [λ 0 ] = M 3 L 3 T −2 , one can check that the above equation (63) has the desired squaredenergy dimension.In order to get an equation involving only dimensionless parameters, one can simply divide by the squared energy associated to the motional quanta (ℏω x ) 2 , which allows us to identify the following dimensionless renormalized parameters as well as the following dimensionless bare couplings Finally, one can rewrite the trapped-ion tadpole selfconsistency Eq. ( 63) in terms of dimensionless quantities as which has the same mathematical structure as the lattice self-consistent equation previously found (58).Again, Matsubara mode sums are to be performed analytically attending to (59).The only difference, apart from the pre-factor (l/d) 3/2 (η N1 (1)) 1/2 , is that the part that depends on the spatial momentum ( k/ω x ) in the propagator (61) contains now the dipolar terms (62).Following this procedure, we can find the trapped-ion analogue of all the selfconsistency equations that include the sunrise diagram (55), which now depend on dimensionless parameters that can be numerically adjusted.
The numerical routines introduced in the previous subsection can be applied to the trapped-ion case directly.We can now use this numerical simulation to approximate the value of the physical mass (55) as the trap frequencies are modified, and the trapped-ion chain gets closer to the linear-to-zigzag phase transition.With this value, we can estimate the effect of the interactions on the effective Compton wavelength (67) ξ eff → ξ eff,P , as well as on the spin-spin coupling strength J eff → J eff,P .The former is due to the quantum and thermal contributions to the physical mass (53), whereas the later comes from the field rescaling with the wavefunction renormalization (54).Both expressions would enter in a renormalised version of the spin-spin couplings of Eq. ( 29), as we recall that these are mediated by the excitations of the selfinteracting scalar field, which are controlled by the physical pole and the full propagator.In particular, we find In this way, one sees that quantum and thermal effects in the λϕ 4 model will change the spin-spin interactions of the trapped-ion quantum simulator.This opens a very interesting perspective, allowing future experiments to probe the nature of the fixed point of this effective QFT by measuring the dynamics of the spins.In fact, the distance dependence of the spinspin couplings has been inferred using various experimental techniques in recent years [43,45,46].Using these techniques while gradually approaching the linear-to-zigzag phase transition would allow one to infer the flow of the critical point, and address universal properties of the QFT in a quantum simulator.For the sake of completeness, we also note that the renormalization of the quartic coupling (20), and the associated four-point functions, would give rise to four-and higher-spin interactions [53].These are, however, much weaker and negligible in a trapped-ion experiment given the constraints considered in this work (32).
To make the numerical results closer to the trappedion language, one can also make use of the Bose-Einstein distribution (33) to obtain a contour plot of the mean number of phonons n(π/d) for the relevant zigzag mode.The final result is that of Fig. 7 where, rather than plotting the Compton wavelength as a function of the dimensionless bare mass, we plot it as a function of the radial trap frequency, which is the standard experimental parameter used to control the shape of the ion crystal.We can see in Fig. 7 (a) how the classical critical point (15) flows with the temperature and with the radial trap frequency (blue dotted line).In the coloured region, which lies well within the symmetry-preserved phase (i.e.linear ion chain), we see how the effective Compton wavelength (67) entering the spin-spin couplings (29) changes as one approaches the critical line.Therefore, using some of the experimental techniques based on probing the real-time dynamics of the spins [43,45,46], one can extract the distance decay and test how the effective Compton wavelength gets renormalised by quantum and thermal effects.In Fig. 7 (b), we also represent the average number of phonons in the zigzag mode as a function of the temperature and the transverse trap frequency.This sets the target detuning of the laser cooling on the ion crystal used to control the contribution of the thermal masses, and the renormalization of the spin-spin couplings in an experiment.
Note that, in order to comply with the constraint (32) and still have spin-spin coupling strengths that are not too slow in comparison to additional experimental sources of noise such as dephasing or motional heating/decoherence, one should avoid getting extremely close to the structural phase transition.There, the renormalised version of the zigzag mode (13), which is proportional to the physical mass m P of the QFT, softens ω zz,P → 0. In light of the constraint in Eq. (32), the laser beatnote ∆ω L < ω zz,P must be very small, leading to very slow spin dynamics.For this reason, our plot in Fig. 7 restricts to the colored regions, and does not consider all the parameters down to the critical line.In fact, the renormalised zigzag frequency, which changes according to Fig. 7(b) is always higher than ω zz,P /2π ≥ 318 kHz, leaving enough frequency space for the spindependent dipole force to fulfill Eq. (32) and still lead to sufficiently-fast spin dynamics in the 0.5-10ms scale.Altogether, Fig. 7 provides a quantitative prediction of how the range of the spin-spin interactions changes as a function of the temperature, and recall that this is a direct consequence of the underlying interactions and Feynman diagrams of the coarsegrained λϕ 4 QFT, and cannot be accounted for if one truncates the description of the system at the quadratic order for the phonons.

Conclusions and outlook
In this manuscript, we have presented a self-consistent approach to estimate non-zero temperature effects in the trapped-ion quantum simulators of spin models.We have argued that the range of the spin-spin interactions mediated by the transverse phonons of the ion chain can be accurately captured by an effective QFT of a real scalar field that has a Yukawa-type coupling to the spins.In the vicinity of a linear-to-zigzag transition of the ion chain, λϕ 4 interactions must be included in this QFT, which can modify the nature of the Yukawa interactions through phonon-phonon scattering.In light of the renormalizability of this QFT, this interaction effects can be recast in a renormalization of the bare quartic coupling, bare mass, and wavefunction renormalization.We have argued that the later two effects yield a renormalization of the range and magnitude of the Yukawa-type spin-spin couplings J ij , which could be inferred from trappedion experiments that reconstruct J ij from the realtime dynamics of the spins.Accordingly, the trappedion quantum simulator could be used to probe renormalization of this paradigmatic QFT.
To find a quantitative prediction of these effects, we have presented a self-consistent approach that resums tadpole-like Feynman diagrams to all orders of the quartic coupling (42).This so-called self-consistent Hartree approximation is afflicted by an infra-red divergence when trying to use it to determine thermal and quantum effects of the critical point of the linearto-zigzag transition.We have thus extended our approach beyond mean-field theory by also considering the sunrise diagram, and its additive and multiplicative contributions to the renormalized mass of the QFT (55).We have discussed how these selfconsistent approach can be applied to the trapped-ion case at non-zero temperatures, which requires using a specific propagator (61) that includes a dipolar regularization (62) of the QFT stemming from a multipole expansion of the Coulomb interactions among the trapped ions.Moreover, we have also discussed how to apply the Matsubara formalism in Euclidean time to account for non-zero temperatures in the experiment, i.e. laser cooling to a non-zero mean number.Using realistic parameters for recent experiments with long ion chains [150], we have been able to derive specific quantitative prediction of thermal effects in the spin-spin interactions, which have been depicted in Fig. 7.As an outlook, we note that these predictions can serve as a guide to future trapped-ion experiments that aim at exploring the present connection between spin-model quantum simulators and this relativistic Yukawa-type problem.We note that the experimental quantum simulation, when working in the long-wavelength regime discussed in this work, will actually go beyond our approximations, and effectively compute the contributions to the Compton wavelength stemming from all possible Feynman diagrams.Moreover, the scaling of this quantity with the experimentally-tunable parameters will unveil the critical exponents of the phase transition of the λϕ 4 model, which should correspond to those of the Ising universality class.Finally, we would like to mention that an interesting problem for future study would be to go beyond the schemes for the spin models of the Ising type, and consider further terms that can lead to an effective relativistic QFT of Dirac fermions Yukawa-coupled to the self-interacting scalar field.This trapped-ion quantum simulator would be closer to a lower-dimensional version of the fermion-Higgs sector of the electroweak interactions, and can provide a way to go beyond semi-classical calculations for the fractionalization of charge in fermion-scalar QFTs [181].

A Thermal effects and Matsubara mode sums
In this Appendix, we evaluate the thermal corrections resulting from the tadpole and sunrise diagrams in the lattice regularized λϕ 4 QFT.The contribution of the resumed tadpole family to the self-energy is The sum over the Matsubara frequencies runs on n 0 ∈ Z and n 1 = 1, .., N 1 , with N 1 the number of sites of the lattice.The lattice analogue of the spatial momentum k is given in (44).The Matsubara sum can be explicitly performed with the help of Cauchy's theorem [159], obtaining where This can be rewritten as Being independent of the external momenta, the tadpole diagrams represent a shift on the bare mass.The first term in the parenthesis is the zero temperature contribution.The second term, which is the mean field analogue of (25), contains the thermal effects.As explained in Sec.3.2, the sunrise diagram contributes both to the mass and wave function renormalization.The finite temperature lattice version of the mass shift (49) is given by T 2 a 2 N 2 1 n0,n1 l0,l1 1 (2πT n 0 ) 2 + ω(n 1 ) 2 1 (2πT (n 0 + l 0 )) 2 + ω(n 1 + l 1 ) 2 1 (2πT l 0 ) 2 + ω(l 1 ) 2 . ( The Matsubara sums can again be performed explic-itly, with the result where we have defined ω 1 = ω(n 1 ), ω 2 = ω(l 1 ) and ω 3 = ω(n 1 + l 1 ) and ω ij = ω k with k ̸ = i, j and i, j Unfortunately, we have not found an analytical expression for the Matsubara sums.They have to be evaluated numerically, which implies truncating the sums to a finite domain n 0 , m 0 = −N 0 , • • • , N 0 .It is thus important to analyze the convergence of the truncated sums with N 0 .For the sake of completeness, we compare in Fig. 8 the numerical evaluation of the tadpole and sunrise mass shifts with the analytical results Eqs.(69) and (73).Whereas the sunrise contribution converges very fast for moderate N 0 at the temperatures of interest, an accurate estimate of the tadpole shift requires much larger values.Therefore, incorporating these expressions in the numerical routines that solve the self-consistent equations makes them much more efficient.The result of the numerical evaluation of the wavefunction renormalization contribution (74) is shown in Fig. 9.As it was the case for the sunrise mass shift, the convergence is very good even in the low-temperature limit for moderate values of N 0 .

B Analytical estimate of the critical ratio
We evaluate here the T = 0 critical point of the λϕ 4 theory in the continuum.As explained in text main text, it is customary to describe the critical point in terms of the dimensionless ratio λ 0 /µ 2 , with µ 2 UVfinite tadpole renormalized mass (45).Within our self-consistent approach, the critical point is determined by Eq. ( 48) and (49) Rescaling the momenta in the integral p, q → µp, µq, we obtain with We evaluate this integral making use of Feynman parameters [148], namely I = d 2 kd 2 q 1 k 2 + 1 (78) Integrating over q we obtain Finally, we use Feynman parameters again to integrate over p, with the result where we have introduced the PolyGamma functions defined in terms of derivatives of Euler's gamma function Γ(z) = ∞ 0 dtt z e −t , namely ψ n (z) = d n+1 log Γ(z)/dz n+1 .

C Critical line crossings
In this Appendix, we describe certain limitations of the current approach that can become important  away for the regime of interest discussed in the main text, namely that of large couplings and masses.As it can be seen in Fig. 10, the critical lines separating the broken and unbroken phases have crossing points in the (m 2 0 a 2 , λ 0 a 2 ) plane.Actually, the critical lines of any two temperatures always cross for sufficiently negative m 2 0 and large λ 0 .This behavior is an artefact that stems from the approximations underlying our procedure, as not only many Feynman diagrams are being discarded in the loop expansion, but also the self-consistency resummations make tadpole-like diagrams prevail upon the rest.Fig. 11 shows that these crossing are however not manifest when the critical lines are plotted as a function of the tadpole renormalized mass µ 2 .In Fig. 12, we also present the contour plots of the physical mass as a function of (µ 2 a 2 , T a), which show a very regular behavior for large values of the couplings and masses.On the other hand, when displaying these contour plots in terms of the bare mass, we encounter an unphysical re-entrance of the symmetry-broken phase at low temperatures, that becomes more evident as one increases the quartic coupling (see Fig. 13).This reentrance is again a consequence of the aforementioned crossings of critical lines.
On the contrary, it is immediate to see that crossings are absent when the critical values of the coupling are plotted as a function of µa.From Eq. (80), such a crossing would imply Γ sr,1 = Γ sr,2 .Since the sunrise mass shift is a monotonic function of the temperature (see Fig. 3), this condition is never satisfied.

Figure 1 :
Figure 1: Trapped-ion chain and zigzag ladder: The equilibrium positions have been calculated numerically solving the non-linear system of equations (2) for a chain with N = 30 ions, to illustrate the linear and zigzag configurations.(a) Linear-chain configuration for a trap-frequency ratio κz = (ωx/ωz) 2 = 10 −4 .(b) Zigzag-ladder distribution for a trap frequency ratio of κz = 0.02.

QFTFigure 2 :
Figure 2: Spin-spin interactions and Yukawa-type QFT: We represent the spin-spin interactions Jij in log-scale for half of a chain of N1 = 30 ions (the remaining half is related by inversion symmetry).The green bars represent the exact expression (31), whereas the red lines stand for the coarsegrained approximation of Yukawa-type interactions mediated by an effective Klein-Gordon field (29).

Figure 3 :
Figure 3: Tadpole and sunrise mass shifts: (a) The tadpole renormalised mass m 2 0 → µ 2 = m 2 0 + Σ td can be expressed in terms of a dimensionless mass shift Σ td /λ0 independent of the coupling.We set µa = 1, and depict both the analytical expression for the T = 0 mass shift (46) (black dotted line), obtained in the thermodynamic limit N1 → ∞, and the finite-temperature expression (red dashed line), obtained by performing the Matsubara sum explicitly (59) for a lattice with N1 = 30 sites.(b) The sunrise renormalised mass µ 2 → mP 2 = µ 2 + Σsr(0) is proportional to a dimensionless ratio µ 2 Σsr(0)/λ 2 0 that is, again, independent of the quartic coupling.We compare the exact T = 0 result in the thermodynamic limit (black dotted line), discussed in Appendix B, with the finite-temperature expression obtained by performing the Matsubara sums explicitly (73) (see Appendix A), for a lattice with N1 = 30 sites (red dashed line).

Figure 4 :
Figure 4: Critical lines for the parity-breaking phase transition: We solve numerically the self-consistency equations (45) and (55) for the critical point m 2P |c = 0. We consider a lattice regularization(43) with N1 = 30 sites and use the analytical Matsubara mode sums discussed in Appendix A along with the Routine R1 described in the text.When the temperature increases, the effect of a thermal mass makes the critical line take lower values for the quartic coupling.For each colored line, the upper region corresponds to the symmetry-preserved ground state (SP-vacuum), whereas the region below corresponds to the ground state with spontaneous breaking of the parity symmetry (SSB-vacuum).

Figure 5 :
Figure 5: Physical mass for the λϕ 4 model on a lattice: Contour plot for non-zero values of the tadpole and sunrise contributions to the physical mass m 2 P a 2 within the symmetric ground state (coloured region).The physical mass is obtained by solving the self-consistent equations (45) and (55).The numerical solution uses the Matsubara mode sums discussed in Appendix A and the Routine R2, described in the text, considering N1 = 30 lattice sites.The dashed arrow indicates the phenomenon of symmetry restoration by increasing the temperature.

Figure 6 :
Figure 6: Critical ratio as a function of λ0: We use Routine R1 and N1 = 2000 to compute λ0a 2 and fc for each µ 2 a 2 value.Three different logarithmic models are fitted (blue line, green dashed line and red dotted line respectively), obtaining f1,c = 20.10(4),f2,c ≃ f3,c = 20.10(5).The inset shows how the the critical ratio f3,c converges to the continuum for different N1 values to check the validity of the discretization used.

Figure 7 :
Figure 7: Renormalization of the range of the spin-spin interactions in a trapped-ion chain: (a) Contour plots of the effective Compton wavelength ξ eff,P /d that controls the range of the Yukawa-type spin-spin interactions.In the expression of the spin-spin couplings Jij(29), one substitutes J eff → J eff,P and ξ eff → ξ eff,P , due to the quantum and thermal contributions in Eq.(67).Thermal effects are encoded in the dependence of the effective Compton wavelength with the temperature.In the x axis, we represent the radial frequency, while the y axis is reserved for the dimensionless temperature(64).We use analytical Matsubara mode sums and N1 = 30 spatial points.Experimental parameters for these simulations have been considered according to a string of40 Ca + ions[150], and fixing the axial trap frequency to ωx/2π = 0.45 MHz, and the laser detuning with respect to the qubit transition to ∆ωL/2π = 0.318MHz.The explored region in (ωz, T ) parameter space avoids getting extremely close to the structrural phase transition, here represented by the blue dotted line, since the constraint(32) for red detunings ∆ω L < ω(k) would imply very slow spin dynamics (see our discussion in the main text).(b) Mean number of phonons n(π/d) = (e (ℏω zz,P /k B T ) − 1) −1 , where ω zz,P refers to the zigzag-mode frequency.(c) Renormalised zigzag mode frequency ω zz,P = mP ωx.

Figure 8 :
Figure 8: Convergence tadpole and sunrise Matsubara mode sums for the mass shifts: Convergence of the Matsubara sums displayed in Fig. 3 as a function the number of modes N0, both for the tadpole (a) and the sunrise contributions (b).

Figure 9 :
Figure 9: Convergence of the Matsubara mode sums for the wavefunction renormalization: Convergence of the Matsubara sums as a function of the number of modes N0 for the quantity ∂ k 0 2 Σsr(0)/λ0 2 appearing in the wavefunction renormalization contribution.

Figure 10 :
Figure 10: Critical lines with crossings: Two critical lines for different temperatures show crossings in the (m 2 0 a 2 , λ0a 2 ) plane as a consequence of the approximations in our procedure.

Figure 11 :
Figure 11: Critical lines with no crossings: (a) We solve the self-consistent equations for the critical point m 2 P |c = 0 and plot λ0a 2 with respect to µ 2 a 2 , which leads to non-crossing critical lines.(b) We also plot the bare mass m 2 0 a 2 as a function of µa 2 to confirm that no crossings occur.

Figure 12 :
Figure 12: Finite-physical mass of the λϕ 4 model on a lattice with respect to µ 2 a 2 : Contour plot for finite values of the tadpole and sunrise contributions to the physical mass m 2 P a 2 , which are obtained by solving the self-consistent equation (55).The numerical solution uses analytical Matsubara mode sums and N1 = 30 spatial points to calculate the tadpole and sunrise diagrams.The white region corresponds to the symmetry-broken phase, and the coloured region to the symmetric one.As λ0a 2 increases, the critical line m 2 P = 0 separating both phases is folded to the right.As shown, no crossings occur when plotting with respect to µ 2 a 2 .

Figure 13 : 2 :
Figure 13: Finite-physical mass of the λϕ 4 model on a lattice with respect to m 2 0 a 2 : Contour plot for finite values of the tadpole and sunrise contributions to the physical mass m 2 P a 2 , which are obtained by solving the self-consistent equation (55).The numerical solution uses analytical Matsubara mode sums and N1 = 30 spatial points to calculate the tadpole and sunrise diagrams.The white region corresponds to the symmetry-broken phase, and the coloured region to the symmetric one.As λ0a 2 increases, the critical line m 2 P = 0 separating both phases is shifted to the left.Crossings can be appreciated for low temperatures.