Second-order phase transitions and divergent linear response in dynamical mean-field theory

Second-order phase transitions appear as a divergence in one of the linear response functions. For a system of correlated electrons, the relevant divergent response can and does involve many-particle observables, most famously the double occupancy. Generally, evaluating the linear response function of many-particle observables requires a many-particle generalization of the Bethe-Salpeter equation. However, here I show that the divergence of linear response functions in dynamical mean-field theory is governed by a two-particle Bethe-Salpeter equation, even for many-particle observables. The reason for this is that the divergence at the second-order phase transition is produced by the self-consistent feedback of the dynamical mean-field.

Second-order phase transitions appear as a divergence in one of the linear response functions.For a system of correlated electrons, the relevant divergent response can and does involve many-particle observables, most famously the double occupancy.Generally, evaluating the linear response function of many-particle observables requires a many-particle generalization of the Bethe-Salpeter equation.However, here I show that the divergence of linear response functions in dynamical mean-field theory is governed by a two-particle Bethe-Salpeter equation, even for many-particle observables.The reason for this is that the divergence at the second-order phase transition is produced by the self-consistent feedback of the dynamical mean-field.
Electronic correlations lead to a plethora of phases, from metal-insulator transitions [1] and magnetism [2,3] to charge-density waves [3][4][5], Wigner crystals [6][7][8], phase separation [9], superconductivity [10,11], orbital [12,13] and bond order [14,15].Many of these phases are already present in variants of the Hubbard model [16,17].Second-order phase transitions between these correlated phases at finite temperature are of special interest, since divergences occur in the correlation and response functions at these points, associated with the vanishing second derivative of the free-energy functional [18,19].According to the fluctuation-dissipation theorem, the relevant correlation functions are manyparticle observables of higher rank than the order parameter itself.For example, for a ferromagnetic or antiferromagnetic transition, the order parameter ⟨n iσ ⟩ is a single-particle operator, while the relevant correlation function ⟨n iσ n jσ ⟩ is a two-particle operator.
However, this analysis of the Bethe-Salpeter equation appears to be limited to the two-particle correlation functions and thus to single-particle order parameters.This excludes the most simple realization of the metalinsulator transition, where the double occupancy D and its response to a change in the Coulomb interaction strength dD/dU are the quantities of interest [35,41,42], so the order parameter is a two-particle operator and the divergent correlation-response function involves fourparticle operators, whereas the compressibility d⟨n⟩/dµ does not diverge at the critical point of the particle-hole symmetric Hubbard model [38,43,44].
More generally, considering the free energy as a function of µ, U and possible other parameters, thermodynamic stability is a condition on eigenvalues of the second derivative matrix of the free energy [41], which can be expressed in terms of mixed response functions like ∂D/∂µ| U .Finally, in multi-orbital systems, higher-order crystal field and magnetic order parameters [45][46][47][48] do not always have a representation as a single-particle observable, which follows from the addition rules for angular momentum in many-electron systems.
Thus, it is relevant to study the response of manyparticle observables in correlated electron systems, especially with an eye on possible divergences.For the particular case of the double occupancy, Kowalski et al. [41] have used the Galitskii-Migdal formula to reduce the problem to single-particle objects, but a more general and systematic approach is clearly beneficial.
Here, I will show that in DMFT the linear response of many-particle correlation functions and especially their divergence is still governed by the usual, two-particle Bethe-Salpeter kernel.In fact, the many-particle order parameter and applied field only show up as "capping stones" at the end points of the two-particle Bethe-Salpeter equation.Thus, they do not generate the divergence at the second-order transition and their role is restricted to determining if the divergence is picked up in a particular response function.The reason for this remarkable simplification, from many-particle to two-particle physics, can be traced back to the particular form of the DMFT equations, where the self-consistent feedback of the dynamical mean-field is responsible for the secondorder phase transition [18,19,35,36].On the other hand, going beyond linear response, the two-particle Bethe-Salpeter is no longer sufficient, as expected.
Consider a general Hubbard model of the form where a, b are sites in a lattice, α and β are orbital labels Here, t aα,bβ and H loc are Hermitian.For a translationally invariant system, t k,αβ denotes the Fourier transform of t aα,bβ to momentum space.The model is considered at a finite temperature T = 1/β, and factors of T are suppressed in the equations for compactness.
In DMFT, this lattice Hamiltonian is replaced by an auxiliary impurity model with the same local Hamiltonian but with a dynamical hybridization function ∆ ν,αβ , where ν is a fermionic Matsubara frequency.This hybridization might be represented as an (infinite) discrete bath to obtain a Hamiltonian formulation of the impurity, or simply as an action in imaginary time.For now, a hybridization of the form ∆ τ −τ ′ ,αβ c † α (τ )c β (τ ′ ) is used, where ∆ ν has been Fourier transformed to imaginary time.The generalization to Nambu space for superconducting phases is discussed at the end.Given a hybridization ∆ ν,αβ , the auxiliary impurity model can be solved numerically [49] and its time-ordered expectation values are denoted by ⟨•⟩.In particular, DMFT works with the imaginary-time single-particle Green's function g αβ (τ ) = c α (τ )c † β and its Fourier transform to Matsubara frequency g ν,αβ .In the following, t k , ∆ ν and g ν are considered as matrices in orbital space, and • −1 denotes the matrix inverse in this space.
The DMFT loop is closed by a prescription to find the hybridization ∆ ν , the dynamical mean-field, which is given by a set of self-consistency conditions, Here, dk = 1/N k k denotes taking the momentum average, i.e., the local part.Equation ( 2) is a coupled set of equations because the solution of the auxiliary impurity model g ν implicitly depends on ∆ ν ′ also for ν ̸ = ν ′ .Linear response of local observables Linear response considers the change of the expectation value of an operator B to a small perturbation H → H − A X of the Hamiltonian, where A is the magnitude of the perturbation and A and X are called conjugate variables.Two examples introduced above are the density of orbital α, nα = c † α c α and the double occupancy on orbital α, Dα = c † α↑ c α↑ c † α↓ c α↓ , which are conjugate to the chemical potential µ and Hubbard interaction U acting on that orbital, respectively.
As in these examples, and in the spirit of dynamical mean-field theory, I focus here on homogeneous local perturbations, i.e., X = sites i Xi [{c † iα , c iβ }], where Xi is a local operator on site i of arbitrary order.Similarly, only site-local observables B are considered.In that case, in DMFT, it makes sense to identify [50] the expectation value of the impurity model as the relevant quantity, i.e., ⟨B⟩ = 1 imp , which can be measured in the impurity solver.For the linear response to the homogeneous field A, the resulting linear change to a local observable B is the same on all sites, i.e., it is a q = 0 response.More generally [51], it is also possible to consider how ⟨B b ⟩ depends on Âa for any pair of sites a, b, and the corresponding q-dependent response function in momentum space.Similarly, since the perturbation is constant in time, the linear response is also assumed to be timeindependent and the response function has ω = 0.The linear response formalism assumes that no spontaneous symmetry breaking in space or time takes place in response to the field, but second-order phase transitions are visible as a divergent linear response.For single-particle operators A X and B, the DMFT linear response is given by the well-known Bethe-Salpeter equation [20].Here, I show that the approach which was previously used to prove the thermodynamic consistency [52] of the DMFT compressibility can also be used to express the linear response of many-particle observerables in simple terms.Derivation For a local (i.e., impurity) expectation value ⟨B⟩, a change in the parameter A of the local Hamiltonian will lead to both direct changes and indirect changes via the DMFT self-consistent field ∆, This requires the calculation of the change of ∆ with respect to A, which can be determined from the fact that the DMFT self-consistency equation has to be satisfied both before and after applying the field.Restating the DMFT self-consistency, Eq. ( 2), in terms of g −1 instead of g will lead to more compact equations in the end [53].
Here, the square brackets denote that the mean-field ∆ depends on A and the inverse of the local Green's function g −1 depends on A both directly and via ∆[A].f is diagonal in ν, so the ν labels are suppressed to keep the notation compact.
As stated before, the objects g −1 , ∆, t k are matrices in orbital space.The derivative of one of these matrices with respect to another matrix is a rank-4 tensor in orbital space.Furthermore, g and ∆ carry a single fermionic frequency, so the derivative ∂g −1 /∂∆ has two fermionic frequencies, i.e., it is a matrix.It will make sense to interpret these rank-4 orbital tensors as matrices (rank-2 tensors) in a space of orbital pairs, keeping the additional matrix structure in frequency space as well.In this pair space, the usual single-frequency rank-2 orbital objects are vectors.For the matrix inverse in this pair space, the notation • ¬1 is used, while • −1 is reserved for the original orbital space.For matrix derivatives, there is the useful identity To satisfy the self-consistency condition after the infinitesimal change in the external field A, Here, g −1 ν = iν − ∆ ν − Σ ν acts as the definition of the impurity self-energy Σ.The relevant partial derivatives of the self-consistency condition are where so-called bubbles of Green's functions are denoted as χ 0 , these are rank-4 tensors in orbital space.They are diagonal in frequency, since f depends on g and ∆ at the same frequency only.Seen as a bubble, both propagators have the same frequency because the ω = 0 response is being considered.In particular, χ 0,lat is the bubble of lattice Green's functions (at q = 0, ω = 0), χ 0,imp is the impurity bubble (also at ω = 0) and χ0 is their difference, the non-local part of the bubble.The only term connecting different Matsubara frequencies, ∂Σ/∂∆ is related to the impurity vertex [19] F , Note that both F and χ 0,imp are rank-4 tensors in orbital space, so they are matrices in pair space and their product is the matrix product in pair space, i.e., another rank-4 tensors in orbital space.Diagrammatically, this corresponds to contracting two legs of both objects.
Inserting these results into Eq.( 6) and using the pairfrequency space notation (i.e., bubbles and vertices are matrices, derivatives with respect to A are vectors) gives Isolating ∂∆/∂A gives with 1 the unit matrix in pair-frequency space.Finally, Here, ∂Σ/∂A is the connected time-ordered correlator T Acc † − ⟨A⟩ T cc † with the fermionic legs amputated [54], while ∂ ⟨B⟩ /∂∆ is the connected time-ordered correlator T Bcc † − ⟨B⟩ T cc † and χ 0,imp ¬1 corresponds to amputating both its fermionic legs.Both depend on a single fermionic frequency (and ω = 0).Finally, ∂ ⟨B⟩ /∂A is the connected time-ordered correlator ⟨T BA⟩ − ⟨B⟩ ⟨A⟩.The ingredients of equation ( 13) are illustrated in Fig. 1, while Fig. 2 contains a diagrammatic representation of Eq. ( 13) itself, where the geometric series ( 1 − χ0 F ) ¬1 has been expanded up to second order in the nonlocal Bethe-Salpeter kernel [19] χ0 F .
Second-order phase transitions Looking at Eq. ( 13), none of the impurity correlation functions can be responsible for the divergence, since the impurity model is a finite system at finite temperature, whose expectation values are smooth functions of the model parameters.Instead, the inversion in Eq. ( 13) is the origin of divergences [19].Exactly at the critical point, one of the eigenvalues of the nonlocal Bethe-Salpeter kernel χ0 F is equal to 1, so the inverse in Eq. ( 13) is divergent, and the

FIG. 2. Diagrammatic representation of the linear response.
The geometric series is shown up to second order, higher-order terms have additional vertices and particle-hole propagators inserted.Lines with arrows represent the non-local propagator G, a pair of these lines represents a nonlocal bubble χ0 .
associated eigenvector V describes the order parameter of the transition.This can be seen as a matrix generalization of the Stoner criterion, where χ0 describes how many electronic fluctuations are available while F is the effective interaction between correlated electrons.In the pair-frequency space view, the second term in Eq. ( 13) is a scalar product of the form vector-matrixvector. The matrix which is inverted in Eq. ( 13) is independent of A and B, i.e., it is always a two-particle kernel, even when A and B are many-particle operators, and a single eigenvector V describes the divergent linear response of any observables with respect to any applied field.The overlap between the eigenvector V and the two capping vertices ∂ ⟨B⟩ /∂∆ and ∂Σ/∂A determines if a particular response function d⟨B⟩/dA is divergent at the critical point.In particular, symmetry can lead to vanishing overlap, thereby avoiding a divergent response.
The Supplemental Material shows the antiferromagnetic transition [55,56] as an example, where the linear response with respect to the staggered field (A = h AF ) is divergent at U = U c while the linear response with respect to the interaction strength (A = U ) is not.The reason is that for the antiferromagnetic transition, the relevant eigenvector V is spin-antisymmetric and thus has non-zero overlap with ∂Σ/∂h which is also spinantisymmetric, but zero overlap with ∂Σ/∂U which is spin-symmetric.
Another example is the metal-insulator transition in the particle-hole symmetric Hubbard model, where at the critical point d⟨D⟩/dU is divergent [41,42] but d⟨n⟩/dµ is not [43] because of frequency symmetry [19,38,44]: the eigenvector V is frequency-antisymmetric while ∂Σ/∂µ and ∂⟨n⟩/∂∆ are frequency-symmetric at particle-hole symmetry [57], so their overlap with V is zero.Below the critical temperature, the resulting hysteresis region has three co-existing solutions (two stable) with different values of ⟨D⟩, but a single value of ⟨n⟩.
Physically, the reason for any divergence in DMFT is a runaway self-consistent response of ∆ to an external perturbation, and the self-consistency equation governing ∆ only involves single-particle operators.In linear response, taking first derivatives thus leads to two-particle correlations only, explaining why the two-particle kernel is sufficient.
On the other hand, the first non-linear response, d 2 ⟨B⟩ /dA 1 dA 2 , requires a three-particle equivalent of the Bethe-Salpeter equation.It enters through the derivative ∂ 2 ∆/∂A 1 ∂A 2 , which can be isolated from a three-particle equivalent of Eq. ( 5).This equation will contain a three-particle impurity vertex ∂ 2 g −1 /∂∆ 2 .Subsequent higher orders require Bethe-Salpeter equations of higher and higher order.Superconductivity Superconductivity can be described in DMFT and its cluster extensions using the Nambu formulation [58][59][60][61], where the dynamical mean-field also has anomalous components of the form ∆ an.(c † c † + cc), which leads to anomalous components in g and in all vertices.To find instabilities towards a superconducting phase, it is necessary to take these anomalous processes into account in the Bethe-Salpeter equation, even in the normal phase, where it corresponds to the particleparticle channel of the nonlocal Bethe-Salpeter equation, see Refs.[62,63] for a recent discussion.For this situation, the present derivation can be generalized by incorporating a Nambu index into the orbital label, which leads to a treatment of the particle-particle and particlehole channels on equal footing.Diagrammatically, propagators and capping vertices with two incoming or outgoing fermions are then allowed.With this generalization, the conclusions about the nature of second-order transitions in DMFT hold, since the necessary ingredient is that the dynamical mean-field couples to precisely two fermionic operators, regardless of their Nambu index.
Locality The approach presented here is restricted to perturbations and operators which are impurity-local and spatially homogeneous, q = 0.The generalization to com-mensurate q ̸ = 0 is straightforward [51,73].On the other hand, an extension to non-local operators, e.g., the linear response to changes in t k or the identification of bond ordering, requires more work.In the same vein, the response to changes in temperature changes the Matsubara frequencies themselves, which requires an extension of the current formalism.
In conclusion, I have shown that the linear response in dynamical mean-field theory is mainly governed by the two-particle Bethe-Salpeter equation, even when manyparticle observables are considered.In fact, the specific form of the applied perturbation and the studied observable only appears as capping vertices at the two ends of the Bethe-Salpeter ladder.This generalizes previous formulas for the (density,double occupancy)-(µ,U ) response matrix [41] to arbitrary local observables and perturbations.The DMFT linear response functions are equivalent to second derivatives of the free energy [18,29,41], so this result shows that any DMFT second-order phase transition or thermodynamic instability must appear in the nonlocal Bethe-Salpeter kernel.Furthermore, the spatial structure of the equation is entirely captured by the nonlocal Bethe-Salpeter kernel, so all divergent response functions have the same correlation length close to the phase transition.

SUPPLEMENTAL MATERIAL CASE STUDY: ANTIFERROMAGNETIC TRANSITION ON THE CUBIC LATTICE
At low temperature and large interaction strength, the Hubbard model has a tendency towards strong magnetic fluctuations.On a bipartite lattice at half-filling, such as the cubic lattice, this takes the form of antiferromagnetic order.In both dynamical and static mean-field theory, the transition is characterized by the number of solutions to the self-consistent equations, as shown in Fig. 3(a).If the staggered field h = 0 then there is always a solution with staggered magnetization m = 0.However, additional solutions with finite |m| for h = 0 appear at U = U c (black dot) and continue to exist for U > U c (thick red line).The solution with m = 0 is stable up to U c and not beyond.For h ̸ = 0 and U > U c , initially there are multiple DMFT solutions but the one with the staggered magnetization aligned along the field is the unique global minimum of the free energy (red shaded area).For larger h, the other solutions eventually disappear (blue line).This phase diagram is similar to mean-field theory for the Heisenberg model.Note that there is a difference with mean-field theory for the Ising model, since there is a continuum of solutions with fixed |m| at h = 0, instead of a pair ± |m|.
On a bipartite lattice, the staggered field and magnetization can be made uniform by defining a new spin variable ς = (−1) |r| σ, where ↑= − ↓, as shown in Fig. 3(b).This mapping leaves the Hubbard interaction invariant, the staggered field changes to −hc † i,ς c i,ς , but the hopping term turns into −tc † i,ς c j,−ς .This transformation makes it possible to do the DMFT with one single-site impurity model in the ordered phase or in the presence of a staggered field.Even though the hopping is no longer spin-diagonal, the hybridization, local Green's function and self-energy are spin-diagonal even in the presence of the staggered field h since any local process involves an even number of "hops" on the bipartite lattice.In the same vein, the system remains at half-filling at µ = U/2 for a bipartite lattice.The mapping to rotated spin variables to make an antiferromagnetic spin pattern uniform can also be generalized to other q [51].

NUMERICS WITHOUT STAGGERED FIELD
The transition has been studied numerically in previous works [55] and there is a publicly available dataset [56,76] for the case h = 0 for the cubic lattice [77], shown in Fig. 6.The U c at which a finite staggered magnetization m appears decreases with the inverse temperature β, we concentrate on the case β = 5, U c ≈ 4. The leading eigenvalue of the dual Bethe-Salpeter kernel approaches unity when U approaches U c from below, signalling the increasing strength of the selfconsistent feedback in DMFT. Figure 3(d) shows the corresponding eigenvector V σσ ′ (ν n ) for U = 3.6 (eigenvalue 0.9) and U = 3.8 (eigenvalue 0.95).The eigenvector lives in the same vector space as single-particle Green's functions, but does not have all of the properties of a Green's function (no causality guarantee, arbitrary normalization and complex phase), since it denotes possible changes in the Green's function.In the present system, the phase can be chosen so that the eigenvector is entirely real.From Fig. 3(d), it is clear that the leading eigenvector for the antiferromagnetic transition is spin-antisymmetric, i.e., V σσ ′ (ν) = σδ σσ ′ V (ν) and frequency-symmetric, V (ν) = V (−ν).The eigenvector Data for h < 0 obtained from h > 0 by symmetry, the circles indicate the h = 0 benchmark [56,76], see also Fig. 6.For Σ, the lowest Matsubara frequency ν0 = πT is shown.Note that at every h, only solutions that are a global minimum of the free are shown.With the additional metastable and unstable solutions, continuous S-shaped curves or loops would form, see Ref. [42].
is non-zero for all Matsubara frequencies.The eigenvector apparently does not change much between U = 3.6 and U = 3.8, even though the eigenvalue gets substantially closer to unity.In addition to the eigenvector with S z symmetry shown here, by SU (2) symmetry there is a second eigenvector of the Bethe-Salpeter with the same eigenvalue, which corresponds to spontaneous order with in-plane magnetization.This eigenvector is of the form and is therefore also spin-antisymmetric.Since our finite field calculations have h along the z-axis, we focus on the V σσ ′ = V ↑↑ − V ↓↓ eigenvector here.

INTERACTION STRENGTH RESPONSE
We can now consider the response along two directions in parameter space, namely changes in U while keeping h = 0 or an applied field h at some constant U , either at U < U c or U > U c .We start with the latter case.Figure 6 shows additional (dynamical) local observables in addition to the magnetization, all at h = 0. Shown are the double occupancy D = n i,↑ n i,↓ , the difference in the real part of the self-energy between majority and minority, Re(Σ ↑ − Σ ↓ ), evaluated at the lowest Matsubara frequency ν 0 = π/β and the average of the imaginary part of the majority and minority self-energy, Im(Σ ↑ + Σ ↓ ), also evaluated at the lowest Matsubara frequency.We observe that all four observables are continuous but not smooth at U c .The solid lines show fits to the final data points in the disordered phase (m = 0), while the dashed lines are fits to the first data points in the disordered phase (|m| > 0).In the bottom panels of Fig. 6, the disordered phase fit (solid line) is subtracted, to more clearly show the non-smooth behavior around U c .Since all shown observables B are continuous functions of U , the linear response with respect to the interaction, dB/dU , is not divergent, but the non-linear response d 2 B/dU 2 can be divergent.
This can be understood based on the analysis in the main text, since the capping vertex ∂Σ/∂U | ∆ and its overlap with the leading eigenvector of the non-local Bethe-Salpeter equation decides if the linear response is divergent [78].Here, approach the critical point from the disordered phase with h = 0 and U < U c , ∂Σ/∂U | ∆ is spin-symmetric, so the overlap with the spin-antisymmetric leading eigenvector is zero.Thus, there is no divergent linear response of any observable with respect to U .

FIELD RESPONSE
To see the diverging linear response, we have to look at the effect of an applied staggered field h, i.e., the vertical direction in Fig. 3(a).In that case, the capping vertex for the applied field, ∂Σ/∂h, has a spin-antisymmetric and frequency-symmetric component already at the Hartree level, so there is a non-zero overlap with the leading FIG. 5. Frequency-dependence of the self-energy in the presence of an applied staggered field h, νn = (2n + 1)πT .DMFT for the Hubbard model on the cubic lattice at β = 5 (Uc ≈ 4).Data for h < 0 obtained from h > 0 by symmetry, the circles indicate the h = 0 benchmark [56,76].
eigenvector.Figure 4 shows finite-field results.Close to the phase transition at U c ≈ 4, the linear response of all shown observables has a divergence, but there are two qualitatively different cases due to the symmetry of the observables.Starting with the h-antisymmetric observables m and Re(Σ ↑ − Σ ↓ ) we find the expected linear response at U < U c , which turns into a discontinuous jump for U > U c .This jump should be understood in the spirit of a Maxwell construction, where an S-shaped curve would appear if additional solutions were also shown.At higher Matsubara frequencies (not shown), Re Σ behaves qualitatively similar, although there are quantitative difference because the DMFT self-energy is a dynamic object.At large frequency, Re Σ and ⟨m⟩ are related by the asymptotic expression lim ν→∞ Re Σ ς (ν) = U n −ς .
For the h-symmetric observables D and Im Σ(ν 0 ), the field-dependence is quadratic and smooth for U < U c .For U > U c , there is a derivative discontinuity at h = 0.If the metastable and unstable solutions would also be drawn, the curves would form a loop, as in Ref. [42].As for the real part, other frequencies of Im Σ(ν n ) are qualitatively similar, as is visible in Fig. 5. Going to high frequency, the asymptotic relation lim ν→∞ Im Σ ς = U 2 iν ⟨n −ς ⟩ (1 − ⟨n −ς ⟩) at half-filling simplifies using ⟨n −ς ⟩ (1 − ⟨n −ς ⟩) = (1 − m 2 )/4, so Im Σ is continuous but not smooth at h = 0 because the same holds for m 2 .
The divergent response of these observables is explained by the capping vertex ∂ ⟨B⟩ /∂∆| h=0 and its spin symmetry.That is, for all four shown observables, their impurity value changes if ∆ changes along the direction given by the eigenvector in Fig 3(d).Here, it is worth pointing out that we are looking at B = Re Σ(ν) and B = Im Σ(ν) one Matsubara frequency at a time, it is possible to get zero response instead of a divergence if one looks at a linear combination like Re(Σ(ν) − Σ(−ν)), because of the usual symmetry of Σ.In terms of the Bethe-Salpeter analysis, this is because the corresponding capping vertex of this linear combination is zero by symmetry.

COMPARISON TO MOTT TRANSITION
For the non-magnetic Mott metal-insulator transition in DMFT [19,20,38,42], the leading eigenvector is spinsymmetric, i.e., V σσ ′ (ν) = δ σσ ′ V (ν) and frequency-antisymmetric, V (−ν) = −V (ν).In that case, the response with respect to changes in the interaction strengh is divergent at the critical point, but the response with respect to the chemical potential is not divergent [43].For the metal-insulator transtion, the lack of divergence in d ⟨n⟩ /dµ is explained by the frequency symmetry of the eigenvector [44], instead of the spin symmetry.
An important difference is that the order parameter of the antiferromagnetic transition corresponds to a clear spontaneous symmetry breaking.As a result, the line h = 0 in the phase diagram is special and everything is (anti)symmetric around this line.The first-order transition takes place exactly at h = 0 for U > U c .For the Mott transition in the U -T plane, there is less symmetry, so the first-order transition line U c (T ) in the hysteresis region is a curve instead of a straight line.
In principle, for the antiferromagnetic transition, there could be observables B for which d⟨B⟩/dh is neither zero nor divergent at the critical point, similar to the finite d( dn dµ )/dU for the metal-insulator transition.This would happen if the capping vertex d ⟨B⟩ /d∆ has non-zero overlap with dΣ/dh, zero overlap with the leading eigenvector of the nonlocal Bethe-Salpeter equation, but non-zero overlap with other eigenvectors with eigenvalue smaller than 1.The first condition requires a spin-antisymmetric component, so then the second requirement means that the overlap with V should vanish due to the frequency structure (instead of the spin structure).This did not happen for any observables studied here, but in principle one can construct by hand linear combination of observables where this happens.In multi-orbital Hubbard models, the larger number of reasonable observables makes this scenario more likely.

SIMULATION DETAILS
The finite-field DMFT simulations shown here were performed using TRIQS [79], the TRIQS CTHYB solver [80].Close to a second order phase transition, the convergence of the DMFT self-consistency cycle slows down and many iterations are needed [19], here approximately 100 iterations were used.The unit of energy t = 1 is used.Finite field calculations were performed for B = 0.003, 0.005, 0.01, 0.02, 0.03, 0.04 and 0.05.The eigenvectors in Fig. 3(d) are calculated without applied field, using w2dynamics [81] for the vertex that enters the Bethe-Salpeter equation and TPRF [82] for the evaluation of the Bethe-Salpeter equation.

FIG. 3 .
FIG. 3. The antiferromagnetic transition on a bipartite lattice, in DMFT.(a) Sketch of the phase diagram at fixed temperature, h is the staggered field and U the Hubbard interaction.(b-c) Antiferromagnetic order before and after the transformation ς = (−1) |r| σ.(d) Leading eigenvector of the dual Bethe-Salpeter kernel for U just below Uc, the corresponding eigenvalue is 0.9 at U = 3.6 and 0.95 at U = 3.8.

FIG. 4 .
FIG.4.The impact of an applied staggered field h.DMFT for the Hubbard model on the cubic lattice at β = 5 (Uc ≈ 4).Data for h < 0 obtained from h > 0 by symmetry, the circles indicate the h = 0 benchmark[56,76], see also Fig.6.For Σ, the lowest Matsubara frequency ν0 = πT is shown.Note that at every h, only solutions that are a global minimum of the free are shown.With the additional metastable and unstable solutions, continuous S-shaped curves or loops would form, see Ref.[42].
which includes spin), t aα,bβ is the hopping and H local is a local Hamiltonian, which is a function of the creation and annihilation operators c † aα , c aβ on that particular lattice site.The local Hamiltonian includes many-particle terms such as the Coulomb interaction1 arXiv:2401.04042v2 [cond-mat.str-el]14 Jun 2024 ( FIG.1.Diagrammatic representation of contributors to the response.The black lines with arrows indicate fermionic propagators.Note that some of the fermionic propagators are amputated and some are not, see main text for definitions.The operators A and B are denoted by small blue and red dots, respectively.