Universal Quantum Criticality from a Gaussian Variational Manifold

Introduction

Introduction. Fermionic Gaussian states form a low-complexity sector of quantum many-body Hilbert space: they are efficiently specified by covariance matrices and underpin protocols for fermionic learning, simulation, and matchgate quantum information [1, 2, 3, 4, 5]. Recent work has quantified fermionic non-Gaussianity as a many-body and computational resource [6, 7], while developing rigorous tests of closeness to the Gaussian manifold and revealing nontrivial structure even within free-fermionic states [8, 9, 10]. At the same time, free and Gaussian fermions can exhibit universal entanglement phenomena, thermalization, long-range propagation, topology, and critical resource signatures [11, 12, 13, 14, 15]. For interacting systems, Gaussian states have long served as variational manifolds [16, 17, 18, 19, 20, 21] and are increasingly used as starting points for more general constructions [22, 23, 24, 25, 26]. Related measurement and decomposition schemes use their exact structure to access interacting fermionic observables [27, 28, 29, 30, 31], while systematic extensions provide controlled routes beyond Gaussianity [32, 33]. These developments raise a sharper question than whether Gaussian states yield accurate energies: must a variational manifold faithfully approximate an interacting critical wave function in order to capture its universal long-distance physics?

The ferromagnetic long-range transverse-field Ising chain, provides a stringent setting for this question. The nearest-neighbor limit is exactly free-fermionic [34, 35], whereas algebraically decaying spin couplings generate nonlocal Jordan–Wigner strings and an interacting fermionic problem. Its ferromagnetic critical behavior spans long-range mean-field scaling, an interacting regime with continuously varying exponents, and eventual recovery of short-range Ising criticality; this structure has been studied by field theory, thermodynamic linked-cluster expansions, tensor-network methods, and quantum Monte Carlo [36, 37, 38, 39, 40, 41, 42, 43, 44]. Entanglement properties of long-range Ising chains have also been studied [45]. Long-range Ising interactions are also directly accessible in programmable trapped-ion quantum simulators [46, 47, 48, 49]. More broadly, Rydberg-atom arrays [50] and ultracold polar molecules [51] provide platforms for engineering long-range interacting spin models.

The precise long-range-to-short-range crossover remains unsettled despite early renormalization-group and large-$N$ analyses [52, 53, 54, 55], subsequent numerical and field-theoretic studies supporting different scenarios [56, 57, 58, 59, 60, 61, 62, 63], and recent large-scale studies of several classical long-range models together with a quantum $O(n)$ renormalization-group analysis [64, 65, 66, 67, 68, 69]. Dangerously irrelevant variables and slow corrections to scaling can complicate finite-size analyses in parts of the phase diagram [43, 44, 70], making a direct thermodynamic formulation especially useful. Previous fermionic-Gaussian work showed that generalized Hartree–Fock states give accurate results for the long-range transverse-field Ising model in the antiferromagnetic regime [21].

Here we develop a thermodynamic Gaussian variational description of the ferromagnetic long-range transverse-field Ising chain and use it to investigate how accurately a noninteracting variational manifold can capture interacting critical behavior. We retain the full spin Hamiltonian and restrict only the variational state to the fermionic Gaussian manifold. Translation invariance and Gaussian purity reduce the variational problem to a momentum-space phase, while the long-range spin interactions remain encoded nonlinearly through determinant and Pfaffian correlation functions. The resulting self-consistent Euler equation is the central object of our analysis. Criticality appears directly as a zero of its Euler field in the thermodynamic limit, yielding an analytic self-consistent condition whose numerical solution closely agrees with independent benchmarks. [39, 40, 41, 43, 44] The nonanalytic momentum structure of the same equation, combined with the asymptotic assumptions specified below, provides closed-form analytic expressions for continuously varying critical properties.

This agreement coexists with a reduction in microscopic wave-function accuracy near the transition. Exact diagonalization shows that the optimized Gaussian state is least faithful to the exact ground state in the critical region, where the exact state also exhibits enhanced non-Gaussianity. Thus, the accuracy of the critical predictions does not track the microscopic state fidelity: the optimized Gaussian description can retain the long-wavelength structure governing critical behavior even where its wave-function approximation is weakest. These results identify the singular structure of the variational solution as the key to understanding how a Gaussian manifold can capture critical properties of an interacting system.

Model and variational formulation. We consider the ferromagnetic long-range transverse-field Ising chain

where $J>0$ sets the interaction scale, $h$ is the transverse field, and $\alpha$ controls the interaction range. The long-range interactions are illustrated in Figure 1(a). For the finite periodic systems used in the numerical calculations, the power law is periodized by summing over all lattice images, so that its fixed-distance thermodynamic limit is $J_r=J/r^\alpha$; the explicit finite-size coupling is given in SM Sec. S5.

The nearest-neighbor limit of Eq. (1) maps under the Jordan–Wigner transformation to a quadratic fermionic Hamiltonian with an exactly Gaussian ground state, whereas the long-range spin interaction generates non-local fermionic strings. Our strategy is not to discard these strings or replace the Hamiltonian by a quadratic one. Instead, we retain the full spin Hamiltonian and restrict only the variational state to a fermionic Gaussian manifold. The nonquadratic character of the spin problem then survives through determinant and Pfaffian expressions for its correlation functions. The finite-size Gaussian construction, our Jordan–Wigner conventions, and the corresponding Pfaffian formulas are given in SM Secs. S3 and S4.

The apparent many-parameter Gaussian problem simplifies strongly in the thermodynamic limit. Using the Jordan–Wigner fermions defined in SM Sec. S4, translation invariance makes the mixed Majorana contraction

a Toeplitz matrix specified by a momentum-space symbol $G(q)$. Purity of the Gaussian state constrains this symbol to unit modulus, so that the variational state can be parametrized by a single real phase $\varphi_q$,

This reduction to one function of momentum is central to the analytic treatment below: the many-body optimization problem becomes a variational problem for the phase profile $\varphi_q$, while the interacting spin structure remains encoded in its real-space determinants.

The Fourier coefficients of the symbol are

They directly determine the spin observables entering the energy. In the convention used here, $\langle\sigma_j^z\rangle=-g_0$, whereas the order-parameter correlator takes the Toeplitz-determinant form

Thus even though the trial state is Gaussian, the long-range interaction enters through determinants whose size grows with the spin separation $r$. This is the source of the nonlinear structure of the resulting Gaussian variational theory.

Combining the transverse magnetization and all magnetic correlators, the thermodynamic variational energy density becomes

All dependence on the variational state is therefore carried by the single phase $\varphi_q$, but through the nonlinear Toeplitz determinants $X_r[\varphi]$. Minimizing this functional leads to the self-consistent Euler equation that forms the starting point of our critical analysis.

Variational Euler equation. The energy functional (6) is nonlinear in $\varphi_q$ through the Toeplitz determinants $X_r[\varphi]$. At a variational minimum, $\delta e/\delta\varphi_q = 0$ for every momentum, which can be written as $A_q\sin\varphi_q - B_q\cos\varphi_q = 0$, with

where $\Delta_{ab}^{(r)} = \partial\det T_r/\partial[T_r]_{ab}$ is the corresponding Toeplitz cofactor. It is useful to combine these two quantities into the complex Euler field

The stationarity condition aligns the phase of the Gaussian symbol with this field. Selecting the branch that minimizes the variational energy gives the self-consistency equation

Equations (7)–(10) are exact within the translation-invariant Gaussian manifold. Their functional derivation is given in the End Matter. Criticality occurs when the Euler field itself vanishes at the symmetry momentum, which provides the starting point for the thermodynamic critical-line condition below.

Thermodynamic critical line. Along the ferromagnetic branch continuously connected to the nearest-neighbor Ising critical point, the critical zero of the Euler field occurs at the Brillouin-zone edge $q=\pi$. At this symmetry momentum $B_\pi=0$ identically, so the criticality condition $\mathcal{F}_E(\pi)=0$ reduces to $A_\pi=0$.

Evaluating the cofactor kernel at $q=\pi$ introduces the alternating projection

where the cofactors $\Delta_{ab}^{(r)}$ are determined by the same Gaussian symbol $G$. The condition $A_\pi=0$ then gives

with $G_c$ the self-consistent critical symbol. Equations (10) and (12) therefore form a coupled nonlinear thermodynamic problem for $G_c(q)$ and $h_c$. The transition is obtained directly from the variational Euler equation, without invoking finite-size crossings.

We solve the coupled equations numerically in momentum space and extrapolate the interaction-range cutoff $r_{\max}\to\infty$; implementation details are given in the End Matter. Exact-diagonalization results for the fidelity at $L=20$, shown in Figure 1(b), reveal that the Gaussian approximation is least accurate near the critical line, where the exact ground state also shows its strongest departure from Gaussianity (see SM). Nevertheless, the resulting thermodynamic critical line, also shown in Figure 1(b), agrees closely with independent finite-size estimates from the entanglement-entropy peak and Binder-cumulant crossings. A comparison with quantum Monte Carlo results [44], presented in the End Matter [Fig. 4(a)], shows similarly close agreement. Further diagnostics and finite-size analyses are given in SM Sec. S5.

Analytic magnetic correlation exponent. The same critical Euler field that determines $h_c$ also controls the long-distance magnetic correlations. At the self-consistent critical solution, define its magnitude $\mathcal{E}_E(q)=|\mathcal{F}_E(q)|$ and characterize the zero at $q=\pi$ by

Here $\xi_E$ is a static momentum-space exponent of the variational stationarity equation. We seek the magnetic correlation exponent defined by $X_r=\langle\sigma_j^x\sigma_{j+r}^x\rangle\sim r^{-x_\sigma}$.

The analytic extraction of $x_\sigma$ uses three asymptotic ingredients, discussed in detail in SM Sec. S1. First, the alternating projection of the inverse Toeplitz matrix grows linearly with $r$ at criticality. This behavior is exact at the short-range critical point and for the canonical single-jump Fisher–Hartwig symbol, and we use its continuation along the optimized critical branch. Together with the exact cofactor identity, this gives

Second, the momentum-resolved cofactor kernel has the scaling form $K_r(k)/K_r\to\Psi(kr)$ in the joint limit $r\to\infty$, $k\to0$ with $kr$ fixed. The canonical single-jump problem admits an exact Cauchy-matrix representation that supports this scaling and fixes a nonzero linear small-$k r$ term; the corresponding derivation and short-range checks are given in SM Sec. S1.

To insert this real-space scaling back into the Euler equation, let $\mathcal{K}_r(k)$ denote the momentum-resolved version of $\mathcal{K}_r$, with $\mathcal{K}_r(0)=\mathcal{K}_r$. Using the critical condition $\mathcal{F}_E(\pi)=0$, the field near the Brillouin-zone edge can be written as

This form separates the critical cancellation from the momentum-dependent correction.

There are two contributions to its small-$k$ behavior. The large-$r$ tail of Eq. (15), together with Eq. (14) and the scaling form above, produces a nonanalytic term, while the regular fixed-$r$ part admits an ordinary Taylor expansion:

The Mellin analysis leading to the last line, including the conditions for the convergence of the long-range contribution, is given in SM Sec. S1. Provided the corresponding amplitudes are nonzero and do not cancel, the smaller power controls the magnitude of the Euler field, so that

Thus the continuously varying branch is controlled by the nonanalytic large-distance tail, whereas the short-range branch is controlled by the regular linear term.

The third ingredient concerns the phase rather than the magnitude of the Euler field. Because of the stationarity condition (10), normalization removes the magnitude of $\mathcal{F}_E(q)$ and transfers its phase directly to the critical Gaussian symbol $G_c(q)$. A zero of $\mathcal{F}_E$ can therefore generate a singular phase structure in the Toeplitz symbol controlling $X_r$. The magnitude exponent $\zeta_E$ alone does not determine this phase jump. We close the asymptotic problem by following the single-jump Fisher–Hartwig branch continuously connected to the short-range Ising solution. In this branch the shifted Toeplitz symbol fixes the Fisher–Hartwig exponent once the endpoint phase is known, while the remaining phase-locking condition relates that endpoint phase to the Euler singularity. This gives

The first relation is therefore the phase closure of the Gaussian asymptotic theory, while the second is the single-jump Fisher–Hartwig exponent of $X_r=(-1)^r\det T_r$. Exact short-range and canonical single-jump checks, together with the precise residual phase assumption, are given in SM Sec. S1.

Combining Eqs. (17) and (18) therefore closes the problem and gives

The crossover occurs when the nonanalytic long-range power reaches the regular linear contribution, yielding $\alpha_\ast=11/4$. At this point the two powers coincide and logarithmic corrections are not excluded.

The distinction between exact Gaussian identities, the asymptotic closures entering the analytic solution, and their interpretation in terms of the physical long-range fixed point is discussed in SM Sec. S2.

The static Euler exponent can also be related to the physical dynamical exponent. At a critical point governed by the long-range interaction, scale invariance gives

Independently, on the long-range branch Eq. (17) gives $\zeta_E=\alpha+x_\sigma-2$, and therefore $z=\zeta_E$. This identification does not require interpreting the Euler field as an excitation spectrum: the two exponents are governed by the same long-range scaling combination. In the short-range Ising regime, $z=1$ and the regular Euler contribution independently gives $\zeta_E=1$. The scaling argument is detailed in the End Matter.

Figure 2(a) compares the analytic Gaussian prediction for $x_\sigma$ with finite-size estimates from the optimized spin states and with available numerical benchmarks. Where ordinary finite-size scaling applies, $x_{\sigma} = 2\beta/\nu$, and the analytic result reproduces both the continuously varying long-range regime and the crossover to the short-range Ising value $\beta/\nu= 1/8$. Figure 2(b) compares the resulting prediction $z(\alpha) = \zeta_E(\alpha)$ with numerical estimates of the dynamical exponent obtained from QMC, showing close agreement across the range of $\alpha$ considered.

Entanglement and logarithmic coefficient. The critical Gaussian symbol also determines the entanglement properties of the variational state. For a spatial block, the von Neumann entropy is obtained exactly from the eigenvalues of its restricted Majorana covariance matrix using the standard Gaussian-state construction [71, 72]; the explicit formula used in our numerics is given in SM Sec. S5. At criticality we characterize its leading growth by

In the genuinely long-range regime, $c_{\mathrm{eff}}$ parametrizes the logarithmic entanglement growth and is not assumed to be the central charge of an underlying conformal field theory.

To determine this coefficient analytically, one additional spectral identification is required. We assume that, after removal of smooth factors irrelevant to the leading logarithm, the restricted covariance operator belongs to the same asymptotic spectral class as the corresponding single-jump critical fermionic problem. The relation between singular fermionic correlation symbols, Toeplitz asymptotics, and logarithmic entanglement scaling is well established [73, 74, 75, 76]. For the corresponding exactly solvable free-fermion interface problem, the entanglement spectrum and logarithmic coefficient are known exactly [77, 78]. We make the additional spectral identification

which gives

where $\operatorname{Li}_2$ is the dilogarithm. On the long-range branch, $\delta= \sqrt{5-\alpha}-1$, so that $s = \sin[\pi(\sqrt{5-\alpha}-1)]$ for $1 < \alpha< 11/4$. For $\alpha> 11/4$, $\delta= 1/2$ and hence $s = 1$, reducing Eq. (22) to $c_{\mathrm{eff}}^{\mathrm{GA}} = 1/2$.

The covariance-matrix entropy itself is exact within the Gaussian variational state, whereas the thermodynamic expression Eq. (22) relies on the additional spectral identification above. Its derivation and precise logical status are given in the End Matter. The independent numerical extraction of $c_{\mathrm{eff}}$ from block and half-chain entropies is described in SM Sec. S5. Figure 3 compares these finite-size estimates with the analytic Gaussian prediction, showing good agreement with the exact-diagonalization (ED) results.

Discussion. Our results show that universal critical information can survive a severe restriction of the variational manifold even when the wave-function description deteriorates near criticality. The self-consistent Euler equation determines the thermodynamic critical line, while its Brillouin-zone edge singularity together with the asymptotic Toeplitz structure determines the magnetic exponent. Combined with long-range critical scaling, the same endpoint exponent gives the dynamical relation $z = \zeta_E$, and, with an additional spectral closure, the logarithmic entanglement coefficient. Exact diagonalization makes the contrast particularly striking: the exact ground state is most strongly non-Gaussian and the optimized Gaussian state least faithful near criticality (see SM), precisely where the critical properties remain accurately reproduced. The relevant information is therefore not encoded in microscopic wave-function fidelity, but in the singular structure retained by the optimized Gaussian solution. More broadly, our results suggest that the reach of a variational theory of criticality is set by two ingredients: whether the manifold can represent the infrared structure controlling a universal quantity, and whether variational optimization selects that structure. Gaussian manifolds may therefore remain effective not only at Gaussian or mean-field fixed points, but also at interacting critical points whose essential infrared data can be encoded through mass zeros, soft modes, or singular covariance structure. Their limitations should instead arise when the universal physics depends on additional infrared parameters, such as a Luttinger parameter or compactification radius, that are not generically encoded by the Gaussian covariance structure. The relevant question is thus not whether the critical state is itself close to Gaussian, but whether the infrared information controlling the observable of interest is representable within the variational manifold.

Acknowledgments

M.A.R. acknowledges partial support from CNPq and FAPERJ (grant number E-26/210.062/2023). We gratefully acknowledge the resources on the LiCCA HPC cluster of the University of Augsburg, co-funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)–Project-ID 499211671. This work was supported by the research funding program “Forschungspotenziale besser nutzen!” of the University of Augsburg.

End Matter

The thermodynamic Euler equation and the critical-line condition are exact within the translation-invariant Gaussian variational manifold. Building on these exact relations, the magnetic exponent and entanglement coefficient require additional asymptotic closures, specified below. We collect the corresponding analytic steps here, while the numerical implementation and finite-size analyses are given in the Supplemental Material.

Thermodynamic critical line. The variational energy depends on the phase $\varphi_q$ only through the Fourier coefficients $g_n$ and the Toeplitz determinants built from them. From Eq. (4),

where $\Delta_{ab}^{(r)}=\partial\det T_r/\partial[T_r]_{ab}$. Thus the variation of each long-range spin correlator is reduced to the variation of the same phase profile that parametrizes the Gaussian state.

Substitution into $\delta\varepsilon[\varphi]=0$ gives

Although this equation has the form of a momentum-space alignment condition, it remains nonlinear because the cofactors $\Delta_{ab}^{(r)}$ are themselves functionals of the symbol $G$. Separating the coefficients of $\sin\varphi_q$ and $\cos\varphi_q$ reproduces Eqs. (7) and (8). Choosing the minimizing stationary branch then gives the self-consistency equation Eq. (10).

Along the ferromagnetic critical branch, the zero of the Euler field occurs at the Brillouin-zone edge $q=\pi$. At this momentum the sine contribution vanishes identically, since $\sin[\pi(b-a+1)]=0$, and therefore $B_\pi=0$. The criticality condition is consequently reduced to the single real equation $A_\pi=0$, with

where $K_r[G]$ is the alternating cofactor projection defined in Eq. (11). Hence $\mathcal{F}_E(\pi)=0$ gives directly

The condition is scalar at $q=\pi$, but it remains self-consistent: the coefficients $\mathcal{K}_r[G_c]$ are evaluated from the same critical symbol $G_c$ obtained from Eq. (10). The pair $(G_c,h_c)$ is therefore determined simultaneously in the thermodynamic limit. As an independent benchmark of the thermodynamic critical line, we compare the self-consistent Gaussian solution with the QMC parametrization [44] in Figure 4(a). The agreement is close over the full range shown, with the largest deviation occurring in the strongly long-range regime. The Gaussian points are obtained directly from the coupled Euler equation and criticality condition.

Brillouin-zone-edge scaling and magnetic exponent

We now justify the asymptotic steps connecting the critical Euler field to the magnetic exponent. The first ingredient concerns the large-$r$ behavior of the critical-field kernel. Whenever $\mathbf{T}_r$ is nonsingular, the cofactor identity gives

At the short-range critical point the alternating projection is exactly $\boldsymbol{\eta}^{\mathsf{T}}[\mathbf{T}^{(0)}_r]^{-1}\boldsymbol{\eta}=-\pi r/2$. Our first asymptotic closure is that its linear growth persists along the long-range critical branch, $\boldsymbol{\eta}^{\mathsf{T}}\mathbf{T}_r^{-1}\boldsymbol{\eta}\sim-\kappa r$ with $0<\kappa<\infty$. Together with $X_r\sim r^{-x_\sigma}$, this gives the scaling $\mathcal{K}_r\sim r^{1-x_\sigma}$ used in Eq. (14).

To determine how this real-space tail enters the Euler field near $q=\pi$, define

By construction, $\mathcal{K}_r(0)=\mathcal{K}_r$. The second line is our asymptotic scaling closure, taken for $r\to\infty$, $k\to0$ with $kr$ fixed. Inserting this form into the exact critical subtraction (15) gives for the large-distance contribution

For $0<\mu<1$, the change of variable $u=|k|r$ converts the asymptotic tail into

up to an overall nonzero amplitude. This establishes the nonanalytic power $|\mathcal{F}^{\mathrm{LR}}_E|\sim|k|^{\alpha+x_\sigma-2}$ quoted in the main text.

The regular contribution has a different origin. For every fixed $r$, $\mathcal{K}_r(k)$ has an ordinary Taylor expansion about $k=0$, whose leading nonvanishing term is generically linear. The full Euler field therefore contains both the analytic $O(k)$ contribution and the nonanalytic large-$r$ contribution above. Their competition gives Eq. (17), $\zeta_E=\min\{\alpha+x_\sigma-2,1\}$.

It remains to connect this magnitude exponent to the magnetic correlator. By Eq. (10), normalization removes $|\mathcal{F}_E|$ and leaves its phase as the phase of the critical symbol $G_c$. Our phase closure assumes that, after removal of smooth nonvanishing factors, this singular part is represented by a single Fisher–Hartwig jump continuously connected to the short-range Ising branch. Its jump strength satisfies $\delta=1-\zeta_E/2$, while the corresponding Toeplitz determinant decays as $r^{-\delta^2}$. This gives precisely Eq. (18).

On the long-range branch, $\zeta_E=\alpha+x_\sigma-2$ together with $x_\sigma=(1-\zeta_E/2)^2$ gives $\delta^2+2\delta+\alpha-4=0$. The physical root, $\delta=\sqrt{5-\alpha}-1$, yields Eq. (19). For comparison, Figure 4(b) also shows the long-range mean-field (LR MF) result $x_\sigma/2=(3-\alpha)/4$ in its regime of validity, $1<\alpha\leq5/3$. The Gaussian prediction closely follows the LR MF behavior throughout this regime, with their relative difference shown in the inset. The long-range power reaches the regular linear contribution at $\zeta_E=1$, for which $x_\sigma=1/4$ and hence $\alpha_\star=11/4$. At this coincidence of powers logarithmic corrections are not excluded.

Entanglement-spectrum closure. The finite-block entropy is obtained exactly from the restricted Gaussian covariance matrix, as described in the Supplemental Material. To obtain its logarithmic coefficient analytically, we make the additional spectral identification stated in the main text: after removal of smooth factors, the singular part of the critical covariance operator is assumed to lie in the same asymptotic spectral class as a single-jump fermionic problem [73, 74, 75, 76]. We identify its Fisher–Hartwig jump $\delta$ with the spectral parameter of the exactly solvable free-fermion interface problem through $s=\sin(\pi\delta)$.

In the interface problem, the reference spectral parameter $\epsilon$ and the entanglement-spectrum parameter $\omega$ are related by

[77, 78]. A mode characterized by $\omega$ contributes $h(\omega)=\log(1+e^{-2\omega})+2\omega/(e^{2\omega}+1)$ to the von Neumann entropy. Using the corresponding asymptotic spectral density gives

following the exact free-fermion result [77, 78]. The factor of two implicit in this normalization accounts for the two ends of the interval.

The integral can be evaluated exactly [77]. Differentiating twice with respect to $s$ gives

whose integration yields

Together with Eq. (34), this reproduces the analytic coefficient in Eq. (22). On the short-range branch $\delta=1/2$, hence $s=1$ and $I(1)=\pi^2/12$, giving $c_{\mathrm{eff}}=1/2$.

Thus the covariance-matrix entropy is exact within the Gaussian state, whereas the closed thermodynamic coefficient requires the spectral identification (33). Once this identification is made, the remaining spectral integral and its dilogarithm evaluation coincide with the exact free-fermion result of Refs. [77, 78].

Analytic support for the asymptotic closures

Here we collect the analytic steps underlying the three asymptotic ingredients used to extract the magnetic correlation exponent in the main text. It is useful to separate identities that are exact within the Gaussian manifold from statements that require an asymptotic continuation along the optimized long-range critical branch. The first closure can be proved exactly for the short-range critical symbol and, more generally, for a canonical single-jump Fisher–Hartwig family. The second admits an explicit Cauchy-matrix representation in the same family and leads directly to the $kr$ scaling and Mellin power counting. The final relation between the magnitude exponent of the Euler field and the Fisher–Hartwig jump is a phase-locking condition; its short-range limit and its winding structure are exact, while the continuation of the endpoint phase along the long-range branch is the remaining nontrivial closure.

Exact cofactor identity and inverse-Toeplitz scaling

Recall

and let

Writing the cofactor matrix as $\Delta^{(r)}_{ab}=\det(T_r)(T_r^{-1})_{ba}$ gives the exact identity

Thus the large-$r$ behavior of $K_r$ is fixed once the alternating projection of the inverse Toeplitz matrix is known. At the nearest-neighbor critical point the Gaussian phase is $\varphi_q^{(0)}=q/2$. Hence

It follows that

Introduce

so that

Therefore

The remaining Cauchy quadratic form can be evaluated exactly. More generally, define

and let $\mathbf{x}=C_\lambda^{-1}\mathbf{1}$. The equations $C_\lambda\mathbf{x}=\mathbf{1}$ read

Introduce the rational function

Since $F(a-\lambda)=1$ for all $a=1,\ldots,r$, its zeros and poles fix it to

At large $z$,

whereas its definition gives

Hence

For $\lambda=1/2$, Eqs. (44) and (S15) give the exact short-range identity

The same calculation extends to the canonical one-jump family

Its Fourier coefficients are

so that

Using Eq. (51),

Thus the linear growth of the inverse-Toeplitz projection is exact for the entire canonical single-jump family, not only at the nearest-neighbor critical point.

For the optimized long-range critical symbol we write the singular part as a single jump multiplied by a smooth nonvanishing factor. The asymptotic continuation used in the main text is that this smooth dressing renormalizes only the coefficient of Eq. (56), not its power:

Combining this with $X_r\sim r^{-x_{\sigma}}$ and the exact identity (S3) yields

Equation (57) is therefore exact in the canonical singular sector and is the first asymptotic continuation used for the fully optimized long-range symbol.

Momentum-resolved cofactor scaling

Define the momentum-resolved cofactor sum

For the canonical symbol (S17), the Cauchy representation gives an exact expression for the normalized kernel:

Thus the joint large-distance/small-momentum limit is reduced to the finite-section asymptotics of an explicit Cauchy inverse.

The $kr$ structure can already be seen directly from Eq. (60). Setting $k=u/r$, the oscillatory factors depend on the scaled coordinates $a/r$ and $b/r$ only. Inserting the exact product form of $C_\lambda^{-1}$ and applying Stirling asymptotics to the resulting products converts the normalized bilinear form into a finite continuum integral. Consequently,

The smooth factor dressing the canonical jump is assumed not to change the scaling variable $u=kr$; this is the second asymptotic continuation used in the main text.

An exact finite-$r$ identity fixes the first derivative at the origin. Let

The Cauchy matrix obeys

Multiplying by $C_\lambda^{-1}$ and using Eq. (51) in the derivative of Eq. (60) gives

Hence, in the scaling limit $k=u/r$,

The leading regular correction is therefore explicitly linear and has nonzero coefficient in the canonical single-jump sector. For $\lambda=1/2$, Eq. (64) reduces to

which is the exact short-range check quoted in the main analysis.

Using $K_r\sim D_Kr^{1-x_\sigma}$, define

The long-range part of the subtracted Euler field then has the asymptotic form

For $0<\mu<1$, setting $u=|k|r$ gives the Mellin form

Here $\Psi_{+}$ denotes the appropriate side of the endpoint. The lower bound $\mu\ge0$ ensures convergence of the large-$u$ tail when the scaling kernel is bounded. From Eq. (65), $1-\Psi(u)=O(u)$, so near the origin the integrand behaves as $u^{-\mu}$ and is integrable precisely for $\mu<1$. Therefore, provided the Mellin amplitude is nonzero,

For $\mu>1$, the first Taylor moment is convergent and an ordinary linear term in $k$ is allowed and dominates over $|k|^\mu$. At $\mu=1$, the borderline harmonic sum gives the usual possibility of a $k\log|k|$ correction. Thus, in the absence of an exact cancellation between the regular and long-range amplitudes,

The nonvanishing of the relevant amplitude and the absence of a fine-tuned cancellation are the associated nondegeneracy conditions.

Fisher–Hartwig jump and the residual phase closure

The final step concerns the phase of the Euler field. Since

the phase of the critical Gaussian symbol is fixed by the phase, rather than the magnitude, of the Euler field. This distinction is important: the exponent $\zeta_E$ defined from $|F_E(\pi-k)|$ does not by itself determine the jump of $G_c$.

The winding structure of the Toeplitz determinant can nevertheless be made explicit. Because $[T_r]_{ab}=g_{b-a+1}$, the Toeplitz symbol entering the determinant is the shifted symbol

For the canonical family (S17),

Thus the single Fisher–Hartwig jump lies in the winding branch

and its Toeplitz determinant has the power

Equations (S39) and (S40) are exact within the canonical single-jump representation. In terms of the endpoint phase $\varphi_{\pi-}=s\pi$,

What remains is to relate $\varphi_{\pi-}$ to the magnitude exponent $\zeta_E$. If the complex endpoint asymptotic of the Euler field is phase locked to

with the opposite side fixed by the reality symmetry, then Eq. (72) gives

Combining this with Eq. (77) yields

At the short-range Ising point, $\zeta_E = 1$ and $\varphi_{\pi^-} = \pi/2$, so Eq. (78) and Eq. (80) are exact and give $\delta= 1/2$ and $x_\sigma= 1/4$. The single-jump winding relation (S41) is also exact throughout the canonical family. The genuinely residual assumption is therefore the phase-locking continuation (S42) for the fully optimized long-range Euler field. It is stronger than the magnitude statement $|F_E| \sim|k|^{\zeta_E}$ and should not be inferred from that magnitude alone.

With this phase closure, the long-range branch satisfies

whose continuously connected solution is

When the long-range power reaches the regular linear contribution, $\zeta_E = 1$, the Fisher–Hartwig closure gives $x_\sigma= 1/4$ and hence

Thus the value 11/4 is not inserted as a separate crossover assumption; it follows internally from the long-range power counting together with the single-jump phase closure. The latter is the only one of the three ingredients for which the present argument remains a phase-locking continuation rather than an exact result for the fully optimized long-range symbol.

Logical status of the analytic critical relations

It is useful to distinguish three logically different levels in the critical analysis: exact statements within the Gaussian variational manifold, asymptotic closures used to obtain closed analytic formulas, and statements about the infrared fixed point of the microscopic interacting spin chain.

Exact Gaussian identities versus asymptotic closures. The variational energy functional, the Pfaffian/Toeplitz expressions for spin correlations, the cofactor identities, and the stationarity equation

are exact within the chosen Gaussian manifold. Likewise, on the continuous ferromagnetic branch connected to the short-range solution, $F_E(\pi) = 0$ is the exact Gaussian criticality condition. These statements do not rely on replacing the interacting spin Hamiltonian by a quadratic Hamiltonian. At the nearest-neighbor endpoint, the Gaussian manifold contains the exact free-fermion ground state [35].

The closed critical exponents require additional asymptotic information. As discussed in Sec. S1, the derivation of

uses the large-$r$ inverse-Toeplitz behavior and the scaling form of the momentum-resolved cofactor kernel. The Fisher–Hartwig result $x_\sigma= \delta^2$ is exact once the critical symbol belongs to the corresponding single-jump class [73, 74, 75, 76], whereas the relation

requires the additional phase-locking continuation discussed there. Accordingly, the closed expression for $x_\sigma^{\mathrm{GA}}$ should be understood as the analytic solution of this Gaussian asymptotic closure.

Long-range Gaussian fixed point. Writing $\sigma= \alpha- 1$, the conventional continuum description of the long-range Ising transition contains the quadratic kernel

Such long-range quantum fixed-point descriptions have been developed in Refs. [36, 37, 38, 66]. If the microscopic chain flows to this Gaussian long-range fixed point, balancing frequency and momentum gives

Power counting gives the quartic coupling the scaling dimension $y_u = 3\sigma/2 - d$, so that it is irrelevant for $\sigma< 2/3$ in one dimension [36, 37, 38, 52, 66]. Related rigorous mean-field and Gaussian results are known for classical long-range Ising models and for high-dimensional quantum Ising systems [79, 80, 81, 82], although these results do not directly prove the infrared behavior of the present one-dimensional power-law quantum chain. These exponents are therefore properties of the Gaussian long-range fixed point, while their application to the microscopic long-range quantum Ising chain assumes that the lattice model flows to this infrared theory. Dangerous irrelevant variables can also complicate ordinary finite-size scaling in this regime [43, 44, 70].

The Gaussian variational result is nevertheless very close to the conventional fixed-point form. Expanding the closed expression at small $\sigma$ gives

so the two agree through first order and differ only from $O(\sigma^2)$ onward.

Euler exponent and physical dynamical exponent. The physical long-range scaling relation

is logically distinct from the Gaussian Euler relation

on the nonanalytic branch. The former follows from the long-range fixed point and its nonanalytic spatial kernel [36, 38, 66], whereas the latter follows from the Gaussian stationarity equation. Consequently, literal equality of the physical $z$ and the Gaussian $\zeta_E$ requires, in addition, that the corresponding magnetic exponents agree. The comparison in Figure 2(b) is therefore best interpreted as a numerical test of how accurately the static Gaussian Euler singularity tracks the physical dynamical exponent, rather than as an identification of $\mathcal{F}_E$ with an excitation spectrum. Existing QMC results provide a direct benchmark for this comparison [43, 44].

Long-range–to–short-range crossover. Within the Gaussian asymptotic closure, the crossover occurs when the nonanalytic Euler contribution reaches the regular linear term, $\alpha+ x_\sigma- 2 = 1$. Combining this with the single-jump phase closure gives

Writing $\sigma_\star= \alpha_\star- 1$, this is $\sigma_\star= 7/4 = 2 - \eta_{\mathrm{SR}}$ with $\eta_{\mathrm{SR}} = 1/4$, i.e. the Sak criterion [53]. This shifted boundary differs from the unshifted $\sigma_\star= 2$ continuation of the original long-range RG picture [52]. Thus the Gaussian asymptotic theory naturally reproduces the Sak boundary; the value $11/4$ is not inserted as a separate assumption. Its physical interpretation should nevertheless be separated from the internal Gaussian derivation, since the precise long-range–to–short-range crossover remains unsettled. Numerical and field-theory studies have supported the Sak scenario [44, 62, 63], while recent work has also explored alternative crossover structures [68, 69].

Entanglement coefficient. A similar distinction applies to the logarithmic entanglement coefficient. The entropy obtained directly from the restricted covariance matrix is exact for the optimized Gaussian state [71, 72]. The closed thermodynamic formula for $c_{\mathrm{eff}}$ additionally requires the spectral-class identification described in the entanglement-spectrum closure. The relation between singular fermionic symbols, Toeplitz asymptotics, and logarithmic entanglement is well established [73, 74, 75, 76]. Once the additional spectral identification is made, the remaining free-fermion spectral and dilogarithm formulas are the known exact results of Refs. [77, 78]; the nontrivial step in the present theory is the identification of the optimized long-range covariance spectrum with that reference spectral class.

The overall conclusion is therefore independent of claiming exactness of every closed analytic expression: the Gaussian variational manifold captures the singular organization of the critical theory much more accurately than the global wave-function fidelity alone would suggest.

Complete Gaussian/Pfaffian variational manifold

We first formulate the variational construction for a general pure Gaussian state and an arbitrary reference Fock configuration. The real, translation-invariant manifold used for the thermodynamic calculations in the main text is obtained as a specialization of this general framework.

Gaussian pairing state and normalization

We consider a fermionic Gaussian state built on an arbitrary fermionic occupation configuration of $L$ spinless modes. Let $c_i^\dagger$ and $c_i$ ($i=1,\ldots,L$) denote the canonical fermionic creation and annihilation operators, satisfying

We choose a reference Fock configuration

and introduce the particle–hole transformed fermionic operators

By construction,

so that $|C\rangle$ plays the role of the fermionic vacuum for the $d_i$ modes. Starting from the reference configuration $|C\rangle$, we introduce the variational Gaussian pairing state

Here $\mathbf{R}$ is an $L\times L$ antisymmetric matrix, which is the natural structure required for fermionic pairing. Equivalently, this state admits the computational-basis expansion

where

and $\mathbf{R}_{I(C)}$ denotes the principal submatrix of $\mathbf{R}$ restricted to the sites in $I(C)$. The fermionic reordering sign is

We use the convention $\operatorname{pf}(\varnothing)=1$. Since the exponential changes the occupation relative to $|C\rangle$ only in pairs, the state $|\mathbf{R},C\rangle$ has definite fermion parity, equal to that of the reference configuration:

In particular, the state has even parity when $N_C$ is even and odd parity when $N_C$ is odd.

The normalization factor is

For real antisymmetric $\mathbf{R}$, one has

This state provides the variational manifold used in the following and is characterized by the antisymmetric pairing matrix $\mathbf{R}$ together with the reference configuration $C$. The vacuum-reference state $|\mathbf{R},0\rangle$ is recovered by choosing $n_i=0$ for all $i$, in which case $\operatorname{sgn}(0,I)=1$.

Basic fermionic correlators

The Gaussian pairing state (S61) is completely characterized by its equal-time two-point functions. We introduce the normal and anomalous correlation matrices of the particle–hole transformed fermions,

which we collect into $\mathbf{C}$ and $\mathbf{F}$. For the state (S61), one may express these correlators directly in terms of the pairing matrix $\mathbf{R}$ by defining

A standard calculation then yields

so that $\mathbf{C}$ is Hermitian, while $\mathbf{F}$ is antisymmetric, as required by fermionic statistics. To connect with Jordan–Wigner strings, it is convenient to work with the Majorana-like combinations

and the corresponding two-point functions

In terms of $\mathbf{C}$ and $\mathbf{F}$, these matrices can be written compactly as

These objects appear naturally in Pfaffian and determinant expressions for spin correlators in § 2.4. In the common case where the variational optimum can be chosen with real $\mathbf{R}$, one has

Therefore, the off-diagonal elements of both $\mathbf{K}$ and $\bar{\mathbf{K}}$ vanish, and the remaining nontrivial information relevant for $\sigma^x$ strings is encoded in $\mathbf{G}$.

The same contractions also determine the Majorana covariance matrix used below for subsystem entanglement. Since $\gamma_j$ is Hermitian whereas $\bar{\gamma}_j$ is anti-Hermitian, we introduce the Hermitian Majorana vector $\mathbf{a}=(\gamma_1,\ldots,\gamma_L,i\bar{\gamma}_1,\ldots,i\bar{\gamma}_L)^T$, whose components satisfy $\{a_\mu,a_\nu\}=2\delta_{\mu\nu}$. The covariance matrix is defined by

or equivalently $\langle a_\mu a_\nu\rangle=\delta_{\mu\nu}-i\Gamma_{\mu\nu}$. In terms of the contractions introduced above,

For real antisymmetric $\mathbf{R}$, where $\mathbf{K}=\mathbf{I}$ and $\bar{\mathbf{K}}=-\mathbf{I}$, this reduces to $\boldsymbol{\Gamma}=\begin{pmatrix}0&\mathbf{G}^{T}\\-\mathbf{G}&0\end{pmatrix}$. For a pure Gaussian state, $\boldsymbol{\Gamma}^{T}=-\boldsymbol{\Gamma}$ and $\boldsymbol{\Gamma}^{2}=-\mathbf{I}_{2L}$.

Jordan–Wigner representation and spin correlators

We derive here the spin correlators required in the Gaussian calculations. After fixing the Jordan–Wigner and particle–hole conventions, we express the observables in terms of the Gaussian two-point matrices $\mathrm{C}$, $\mathrm{F}$, $\mathrm{G}$, $\mathrm{K}$, and $\bar{\mathrm{K}}$. For real antisymmetric $\mathrm{R}$, the relevant Pfaffian expressions reduce to determinants.

We map spins to the physical fermions according to

so that

The physical fermions and the particle–hole transformed fermions are related by

Introducing $s_l=1-2n_l$, one obtains

and therefore

To keep track of the signs associated with the reference configuration, we introduce the diagonal matrix

For later use, for $l<m$ we also define the reference-configuration string factor

Correlators from fermionic two-point matrices

Here, all spin observables are evaluated from the Gaussian two-point functions. Longitudinal correlators follow directly from $\mathrm{C}$ and $\mathrm{F}$, whereas transverse Jordan–Wigner strings are expressed through the Majorana-bilinear contractions $\mathrm{G}$, $\mathrm{K}$, and $\bar{\mathrm{K}}$ defined in Sec. S3 S3.B. For real $\mathrm{R}$, $\mathrm{K}=\mathbf{I}$ and $\bar{\mathrm{K}}=-\mathbf{I}$, and the relevant Pfaffians reduce to determinants involving only $\mathrm{G}$. A compact summary of the resulting one-, two-, three-, and four-body correlator formulas is given in Table S1.

Correlator	Conditions	Real-$R$ expression
$\langle\sigma_j^x\rangle$	definite parity	$0$
$\langle\sigma_j^z\rangle$	any $j$	$s_j(2C_{jj}-1)$
$\langle\sigma_l^z\sigma_m^z\rangle$	any $l,m$	$s_ls_m[4C_{ll}C_{mm}+4C_{lm}\delta_{lm}-4C_{lm}C_{ml}-4F_{lm}F_{ml}^{*}-2C_{ll}-2C_{mm}+1]$
$\langle\sigma_l^x\sigma_m^x\rangle$	$l<m$	$\chi_{lm}\det\widetilde{\Gamma}^{xx}(l,m)$
$\langle\sigma_l^x\sigma_m^x\sigma_n^z\rangle$	$n<l$ or $n>m$	$\chi_{lm}s_n\det\widetilde{\Gamma}^{xxz}(l,m;n)$
$\langle\sigma_l^x\sigma_m^z\sigma_n^x\rangle$	$l<n<m$	$-\chi_{lm}s_n\det\left[\widetilde{\Gamma}^{xx}_{n-l+1,n-l}(l,m)\right]$
$\langle\sigma_l^x\sigma_m^x\sigma_n^x\sigma_j^x\rangle$	$l<m<n<j$	$\chi_{lm}\chi_{nj}\det\widetilde{\Gamma}^{xxxx}(l,m;n,j)$

One-body longitudinal correlator

Because $\lvert\mathrm{R},\mathrm{C}\rangle$ has definite fermion parity,

Using Eq. (122),

Two-body longitudinal correlator

For $\langle\sigma_l^z\sigma_m^z\rangle_{\mathbf{R},C}$, Wick’s theorem gives

To write the sum over distinct pairs compactly, define

Summing Eq. (127) over $l<m$ yields

For the vacuum reference, $S_C=I$ and $t_C=L$, giving

Two-body transverse correlator

The Jordan–Wigner string in $\sigma_l^x\sigma_m^x$ reduces to the Pfaffian of a $2(m-l)\times2(m-l)$ antisymmetric matrix. Define $\widetilde{\Gamma}^{xx}(l,m)$ through its upper-triangular entries, for $k<k'$,

The remaining entries follow from antisymmetry.

The reference-configuration signs can be collected into

Using $\operatorname{pf}(\mathbf{D}\mathbf{A}\mathbf{D}^{\mathsf{T}})=\det(\mathbf{D})\operatorname{pf}(\mathbf{A})$, one finds

For real $\mathbf{R}$, $\mathbf{K}=\mathbf{I}$ and $\overline{\mathbf{K}}=-\mathbf{I}$, and the Pfaffian reduces to

where

Three-body mixed correlator

For distinct $l,m,n$ with $l<m$, introduce

where

and

The general Gaussian expression is

For real $\mathbf{R}$, the determinant representation depends on the location of the $z$ insertion. If $n<l$ or $n>m$, define

and

Then

If $l<n<m$, the insertion lies inside the Jordan–Wigner interval and

where the indicated row and column are removed from $\widetilde{\Gamma}^{xx}(l,m)$.

Four-body transverse correlator

For $l < m < n < j$, define

where

The correlator is

For real $\mathbf{R}$, define the two cross-correlation blocks

Then

and

Finite-size Gaussian optimization and numerical critical analysis

The thermodynamic Gaussian equations developed in the main Letter provide direct predictions for the critical line, the magnetic correlation exponent, and the logarithmic entanglement coefficient. We now test these predictions independently by optimizing finite-size Gaussian states for the original long-range transverse-field Ising Hamiltonian.

For each system size $L$, interaction exponent $\alpha$, and transverse field $h$, we minimize the expectation value of the finite-size Hamiltonian. Using the spin correlators derived in § 2.4, the variational energy is evaluated directly as

where the power-law interaction is periodized by summing over all lattice images,

with $r_{lm}=\min(\lvert l-m\rvert,L-\lvert l-m\rvert)$. At fixed separation, this prescription has the thermodynamic limit $J_{lm}\to J/r_{lm}^{\alpha}$ as $L\to\infty$. The optimized state is

The finite-size calculation serves two complementary purposes. First, for sizes accessible to exact diagonalization, the optimized Gaussian state can be compared directly with the exact ground state through its energy, fidelity, and fermionic antiflatness. Second, substantially larger systems can be studied through Gaussian correlation functions, allowing independent finite-size determinations of the critical field, magnetic exponent, and logarithmic entanglement coefficient. We expect these quantities to approach, with increasing $L$, the critical behavior obtained independently from the thermodynamic Euler equation in the main Letter.

Variational accuracy and non-Gaussianity

For system sizes accessible to exact diagonalization, we compare the optimized Gaussian state with the exact ground state $\lvert\psi_{\mathrm{ex}}\rangle$. The first diagnostic is the relative energy error, which we define as

Because the variational state is obtained by minimizing the energy, $E_{\mathrm{var}}$ is an upper bound to the exact ground-state energy.

A more direct wave-function comparison is provided by the fidelity

with $0\leq F\leq1$. The fidelity measures how closely the optimized Gaussian wave function reproduces the exact many-body ground state.

To quantify the non-Gaussianity of the exact state itself, we use the fermionic antiflatness introduced in [6]. The Majorana covariance matrix $\Gamma$ is defined earlier in Sec. S3; we therefore do not repeat its construction here. The $k$th antiflatness is

For a pure Gaussian state, $\Gamma^{\mathsf{T}}\Gamma=\mathbb{I}_{2L}$, so that $\mathcal{F}_k=0$. A nonzero antiflatness therefore measures the departure of the target state from the pure Gaussian manifold.

The energy error and fidelity shown in Figure 1 of the main Letter indicate that the Gaussian approximation is least accurate in the vicinity of the transition. Additional finite-size results for the energy error and fidelity are shown in Figure S1, confirming the same behavior across the sizes considered. To probe the departure from Gaussianity directly at the level of the exact many-body state, we additionally examine the fermionic antiflatness in the left panel of Figure S2. Its enhancement near the transition provides finite-size evidence that the exact ground state is most strongly non-Gaussian in the same region where the Gaussian fidelity is most strongly reduced. Thus the accuracy of the Gaussian variational predictions does not rely on the exact critical state becoming Gaussian.

Entanglement peak and finite-size critical field

The half-chain von Neumann entropy provides a first independent finite-size estimate of the critical field. For each optimized Gaussian state, we restrict the Majorana covariance matrix $\Gamma$ to a block of $\ell= L/2$ consecutive sites. Denoting the eigenvalues of the resulting $2\ell\times2\ell$ matrix $\Gamma_\ell$ by $\{\pm i\nu_j\}_{j=1}^{\ell}$, with $0 \le\nu_j \le1$, the von Neumann entropy is

We denote the resulting half-chain entropy by $S_1(L/2;h)$.

For fixed $(L,\alpha)$, the half-chain entropy develops a pronounced enhancement in the vicinity of the transition. We define the location of its maximum by

The finite-size drift of this maximum is then extrapolated to obtain

These extrapolated entanglement-peak positions are one of the two independent finite-size estimates of the critical line used in Figure 2 of the main Letter.

Representative entropy scans are shown in Figure S3. A common feature is the enhancement of $S_1(L/2)$ close to the critical field, followed by reduced entanglement as the system moves farther into either phase. The detailed finite-size shape depends appreciably on the range of the interaction.

For $\alpha=4$, the entropy develops a relatively sharp cusp-like maximum close to the critical field, whose height increases with system size. At $\alpha=2.5$, the peak remains clearly visible but is broader and more asymmetric at finite $L$. For the more strongly long-range case $\alpha=1.6$, finite-size effects are substantially stronger: the entropy maximum is broad, its position drifts visibly with $L$, and the descent on the paramagnetic side sharpens only gradually with increasing size. Consequently, extrapolation of $h_{\mathrm{peak}}(L,\alpha)$ is particularly important in the long-range regime.

These entropy scans therefore play two distinct roles in the analysis: their peak positions provide the finite-size critical-field estimates used in the main critical-line figure, while their values at the extrapolated critical point are used below to determine the logarithmic entanglement coefficient.

Binder-cumulant estimate of the critical field

A second and independent finite-size determination of the transition is obtained from the Binder cumulant. We define the extensive and intensive magnetizations along the ordering direction by

and the Binder cumulant by

The second and fourth moments entering Eq. (161) are computed using the two- and four-spin correlation functions derived in § 2.4.

Because $U_L$ is dimensionless, its critical behavior is less sensitive to nonuniversal amplitudes than that of the order parameter itself. Near the transition, curves for different system sizes become weakly size dependent and exhibit crossings in the vicinity of $h_c$. From the crossing of successive sizes we obtain a size-dependent estimate $h_{\mathrm{cross}}(L,\alpha)$ and extrapolate its drift according to

The values of $h_c^{(\mathrm{B})}(\alpha)$ constitute the second finite-size data set compared with the direct thermodynamic Gaussian critical line in Figure 2. Representative Binder-cumulant curves are shown in Figure S4.

Magnetic finite-size scaling and exponent extraction

The magnetic critical exponent is obtained independently from the finite-size behavior of the squared order parameter

For the long-range model it is useful to write the finite-size scaling form in terms of the generalized scaling dimensions

The scaling forms are

When ordinary finite-size scaling applies, $\theta= 1$, so that $x_m = \beta/\nu$ and $y_h = 1/\nu$. Operationally, the Binder crossings provide an initial estimate of $h_c$. The exponent $y_h$ is constrained by the critical scaling region, while $x_m$ is determined from the collapse and from the size dependence of $\langle m_x^2\rangle$ at criticality.

For a trial set $(h_c,y_h,x_m)$, the rescaled variables are

Representative collapses are displayed in Figure S4. For the values shown there, $\alpha= 2.5$, $3$, and $4$, ordinary finite-size scaling applies and hence $x_m = \beta/\nu$. At $\alpha= 4$ we obtain $\beta/\nu\simeq0.125$ and $h_c \simeq1.156$. For $\alpha= 3$ we similarly find $\beta/\nu\simeq0.125$ and $h_c \simeq1.420$. At $\alpha= 2.5$, the critical field shifts to $h_c \simeq1.759$ and the magnetic exponent to approximately $\beta/\nu\simeq0.1667$, displaying the continuously varying long-range behavior.

The finite-size magnetic exponents obtained in this way provide the numerical Gaussian points compared with the analytic prediction in Figure 2(a) of the main Letter.

Data-collapse estimator. A practical issue in the scaling analysis is the choice of the field interval $[h_{\min},h_{\max}]$. If the interval is too narrow, the fit becomes noisy and unstable, while an excessively wide interval introduces noncritical corrections. We therefore repeat the analysis for several candidate windows.

For each trial window and each set $(h_c,y_h,x_m)$, the rescaled curves $y(x)$ are interpolated onto a common grid of $N_x$ points in their overlap region. Denoting by $y_L(x_k)$ the interpolated curve for size $L$, we define

and use the normalized collapse cost

The optimal parameters minimize $C$, and the preferred field window is the one producing the smallest stable collapse cost.

$R^2$ estimator. As an independent estimate, we analyze the size dependence of $\langle m_x^2\rangle$ at fixed $h$. For every field value we fit

or equivalently

Writing $y_i=\log\langle m_x^2\rangle_i$ and $\widehat{y}_i=a(h)+s(h)\log L_i$, the quality of the fit is measured by

The critical field is identified close to the maximum of $R^2(h)$, and the corresponding magnetic scaling dimension follows from

Thus the preferred fitting interval for this method is the one producing the largest stable value of $R^2$. A representative example of the collapse-cost and $R^2$ diagnostics is shown in Figure S5.

Entanglement scaling and effective central charge

We finally determine numerically the logarithmic entanglement coefficient $c_{\mathrm{eff}}(\alpha)$ compared with the analytic Gaussian prediction in Figure 3. This numerical extraction uses only the finite-size optimized Gaussian states and is independent of the dilogarithm expression derived in the End Matter.

For a subsystem $A$, the covariance matrix is restricted to the $2|A|$ Majorana modes contained in $A$. If its eigenvalues are $\pm i\nu_r$, with $0 \le\nu_r \le1$, we define

The Gaussian Rényi entropy is

while its $q \to1$ limit gives the von Neumann entropy used throughout the main Letter [Eq. (S121)].

For a periodic critical chain described by a local conformal field theory, the block entropy obeys

where $d(\ell)$ is the chord distance and $B$ is nonuniversal. We use this form as a numerical characterization of the logarithmic entropy growth.

There are two complementary finite-size extractions. First, for fixed $L$ we fit the block-size dependence to

which gives

This provides a useful check of the logarithmic form and of finite-size effects.

The main numerical values used in Figure 3, however, are obtained from the half-chain entropy. Setting $\ell= L/2$ in Eq. (176) gives

where the additive term $-(c/3)\log\pi$ has been absorbed into $B'$. We therefore fit the optimized critical Gaussian data to

The entropy entering this fit is evaluated at the corresponding critical field $h_c(\alpha)$. Representative half-chain fits will be shown in Figure S6.

As an additional consistency check, the Rényi entropies may be analyzed using the local-CFT form

If $S_q(\ell) = A_q\log d(\ell) + B_q$, this gives

In a conventional conformal regime these estimates should be independent of $q$ and agree with the von Neumann result.

For sufficiently short-range interactions the extracted coefficient approaches the Ising value $c = 1/2$. For smaller $\alpha$ it varies with the interaction range. In this genuinely long-range regime we do not assume that the coefficient is a Virasoro central charge; we therefore use $c_{\mathrm{eff}}(\alpha)$ as the effective coefficient of logarithmic entanglement growth. The half-chain values obtained from Eq. (180) are the numerical points compared with the analytic Gaussian prediction in Figure 3.