Papers

We investigate Jordan structure and dynamical limits in six families of integrable spin-$\frac{1}{2}$ chains obtained from the boost automorphism method. Local operator algebras connect nilpotent evolution on finite periodic chains to central defects in spectral additivity, polynomial unipotent $R$-matrices and constant braid limits. For conservative endpoint erosion, we determine the characteristic and minimal polynomials and show that the Jordan partition can change while both polynomials remain unchanged. We also identify the deletion limits that distinguish the parity obstruction to removing signs by diagonal similarity in Class 5 from the parity nilpotent dependence in Class 6. Matching expansions and fragmentation recursions yield exact amplitudes and survival probabilities for processes with killing. In addition, we establish criteria excluding conservative Markov realizations and determine which conservative completions of specified reaction densities admit regular $R$-matrices of difference form.

We examine the Casimir interaction between planar massive Dirac materials coupled to a Lorentz-violating Bumblebee-CFJ vector sector. The standard renormalized one-loop polarization tensor of the Dirac quasiparticles is used as the nonlocal surface response in the Lifshitz scattering formalism. For the normally oriented Lorentz-violating background, the surface response remains parity-even. Hence, the repulsive regimes do not originate from a Hall-like conductivity. Instead, they arise from the CFJ-induced mixing of TE and TM bulk modes, combined with the spatially dispersive response of the Dirac sheets. Furthermore, we show that the Dirac mass suppresses the interaction and changes its large-distance scaling, whereas the CFJ parameter drives transitions between attractive and repulsive regimes.

We propose a BRST-invariant formulation of the Gribov–Zwanziger model in the Landau gauge, in which the Gribov scale is produced by the spontaneous breaking of a scalar quartet symmetry, rather than through the explicit breaking term of the conventional construction. We show that the resulting action preserves the ghost equation of the Landau gauge and reproduces the characteristic Gribov-type gauge propagator of the original Gribov scenario, a BRST-invariant mass term of the refined type being accommodated by the same mechanism. We establish the renormalizability of the model to all orders by means of algebraic renormalization techniques, no counterterm being required outside those already present in the classical action. We further show that the condensation of the quartet leaves unchanged the cohomology through which the observables are characterized, the operator realized in the displaced vacuum having the same cohomology as that of the symmetric phase, so that the mass scale is produced at no cost in the physical content of the theory, and the construction of genuine gauge-invariant observables in the Gribov–Zwanziger setting is thereby left open to the usual cohomological methods.

We determine the quadrupolar logarithmic tidal response of electrically charged Einstein-Euler-Heisenberg black holes from the coupled static odd-parity gravitational and electromagnetic equations. Although the second electromagnetic invariant vanishes on the electric background, its quadratic interaction contributes to the axial fluctuation equations. At first order in the Euler-Heisenberg coupling $a$, a direct asymptotic recursion gives the canonical running matrix $C_{\rm can}^{\rm EH}=(9aQ^2/10)\left(\begin{smallmatrix}0&-Q\\-Q&M\end{smallmatrix}\right)$, defined as the coefficients multiplying the canonical source vector in the $r^{-2}\log(r_0/r)$ term. The resulting local coefficients retain the full dependence on the subextremal ratio $\frac{Q}{M}$ at this perturbative order. The two quartic invariants partially cancel in both nonzero channels, while the action-normalized fields make gravitoelectromagnetic reciprocity explicit. The local extremal limit of the running matrix reproduces the corresponding two-invariant effective-field-theory result after convention conversion. Allowing independent quartic couplings further shows that identical electric background geometries can have different axial running matrices and yields algebraic inverse relations within the specified two-operator model. The calculation determines local logarithmic running within minimally coupled quartic electrodynamics; the finite horizon-matched response requires an additional global matching calculation.

In this note we present a direct relation between the geometrical definition of a quasi-Einstein metric and Horava gravity: the solutions of the large-distance effective action of Horava gravity that are static and orthogonal are quasi-Einstein metrics. The converse is also true if the function over the quasi-Einstein manifold satisfies a given equation.

We compute the one-loop breaking of the Slavnov-Taylor (ST) identities in a massive Abelian chiral gauge theory with spontaneous symmetry breaking, regularized in the Breitenlohner-Maison-`t Hooft-Veltman (BMHV) scheme. We also provide the one-loop finite counter-terms required to restore the one-loop ST identities. The recent formulation of spontaneously broken models based on the Fröhlich-Morchio-Strocchi (FMS) gauge-invariant and dynamical (i.e. propagating inside loops) fields is used. The proposed algebraic approach holds true to all orders in the loop expansion and allows for a consistent decomposition of the ST breaking into independent sectors. By using cohomological techniques we eventually constrain the dependence of the ST breaking on the vacuum expectation value of the Higgs field. This yields a set of differential equations that provide a useful check of the mass-dependent terms in the ST breaking terms in the broken phase. They can be generalized to the non-Abelian case, both in the gauge-invariant formalism and in the conventional one. In the massless limit we find agreement with previous calculations present in the literature.

We determine the intrinsic mirror of an orbifold cubic surface with three nodes and interpret the resulting geometry through mirror symmetry between the $PGL(2,\mathbb{C})$ and $SL(2,\mathbb{C})$ character varieties of the once-punctured torus. We further construct a quantization of the mirror using higher genus Gromov-Witten theory and show that it recovers the $A_1$ spherical double affine Hecke algebra.

In a low-energy effective field theory (EFT), scaling logarithmic terms arise from the renormalization of higher-dimensional couplings that absorb divergences generated by loop integrals with lower-order vertices. Consequently, leading logarithms (LLs) are generated entirely by the leading-order Lagrangian, rendering them universal low-energy predictions. Adapting techniques developed for chiral perturbation theory, we analyze the series of LLs in the EFT of general relativity. As a concrete application, we compute LL corrections to the scalar pole mass up to three-loop order and demonstrate that they can be determined recursively to arbitrary loop order. To confirm the physical validity of our results, we carry out the calculation using a general gauge-fixing parameter and explicitly verify its complete cancellation in the physical observable.

The massive Klein-Gordon equation on the Schwarzschild exterior reduces to a radial differential equation of confluent-Heun type. We present its exact reduction by rigorously retaining a subleading centrifugal term that preserves the horizon indicial exponents but critically shifts the accessory parameters. Within the Svartholm–Schmidt hypergeometric-expansion framework, we derive the three-term recurrence relation and establish the corrected continued-fraction compatibility condition matching the minimal tail to the lower-end recurrence. Crucially, we demonstrate that the physical compactification coordinate $z=(r-1)/r$ shifts the irregular singular point to $z=1$, inherently transforming the hypergeometric expansion into a five-term recurrence that cannot degenerate to three terms. To circumvent this obstruction, we construct the horizon-normalized physical branch via a direct Frobenius series at $z=0$. This representation exhibits high-order convergence and near-precision-floor residuals against the exact radial ODE, while clarifying why quasinormal-mode spectral selection remains a connection problem at the irregular singular endpoint.

We study D2-brane holography within the purely dilatonic sector of four-dimensional maximal ISO(7)-gauged supergravity, obtained as a consistent truncation of type IIA supergravity on S$^6$. In the dual frame, the four-dimensional Lagrangian resulting from the reduction of the massless IIA theory exhibits scaling similarity, thus reflecting the generalised conformal symmetry of the dual three-dimensional maximally supersymmetric Yang-Mills theory. Using an auxiliary spacetime construction, we show that the holographic CFT$_3$ fixed points and RG-flows known in the massive IIA theory have natural counterparts in the massless IIA theory: QFT$_3$s with generalised conformal symmetry connected by generalised RG-flows, which we construct numerically. We also revisit the Coulomb branch solutions and propose an extension of Gubser's singularity criterion to supergravity theories exhibiting scaling similarity.

Energy transport across conformal interfaces is universal in two spacetime dimensions, but its higher-dimensional counterpart has remained largely unexplored. In higher dimensions, scattering can depend on additional kinematic data, such as the angle of incidence, and conformal symmetry is far less restrictive. We use gravitational holography to study this problem. For a broad class of AdS$_d$-sliced domain-wall geometries, we derive a closed-form expression for the energy transmission coefficient in terms of the bulk warp factor. Our formula reproduces known two-dimensional holographic results and extends them to arbitrary dimension. Within this class, we establish different aspects of higher-dimensional universality. In particular, we show that the transmission is independent of the incidence angle and of the profile of the incident perturbation. We also derive bounds on transmission controlled by central-charge-like quantities characterizing the number of degrees of freedom of the theories on the two sides of the interface. This provides the first explicit result for reflection and transmission across conformal interfaces in higher dimensions, where no CFT-based results are currently available.

We study holographic spread complexity for a Type IIB solution dual to a deformation of $\mathcal{N} = 4$ SYM that flows to a confining (2+1)-dimensional supersymmetric theory in the IR. We use the proposal that relates the rate of change of spread complexity of local operators with the proper momentum of an infalling probe particle in the bulk. We analyse two cases: a probe particle falling radially and a probe falling radially while rotating with conserved angular momentum. We elaborate on how to define the proper coordinate in the presence of a conserved Noether charge using the Routhian. The trajectories for the probes were studied numerically, revealing an oscillatory behaviour of the complexity.

We consider general EMD-like theories involving arbitrary numbers of both Maxwell fields and dilatonic scalars. We study charged black holes and derive several algebraic and differential relations among the thermodynamic quantities, scalar charges and the non-extremal parameter. These allow us to obtain two sets of master equations of thermodynamic quantities, depending on the choice of basic variables parameterizing the system. Solutions to either set of the master equations give rise to all the thermodynamic quantities, without constructing black hole solutions. We obtain all the thermodynamic quantities perturbatively in the power series expansions of charges, which we then verify using perturbative black hole solutions up to and including the quartic order in charges.

We develop a bootstrap approach for weakly broken gauge theories in (anti-)de Sitter space, focusing on cases in which the conservation of the dual boundary current is broken by a double-trace operator. Integrating the corresponding Ward identity, we derive pseudo-charge conservation identities that contain fixed nonlocal contributions as integrals of three-point functions. These identities impose consistency conditions on the space of allowed theories. As an illustrative example, we study Yang-Mills theory in AdS with symmetry-breaking boundary conditions for charged matter and show that the pseudo-charge conservation identities constrain the mixing between the two quantization sectors. We then apply the framework to the Higgs mechanism for gravity in two AdS spaces with a common boundary, coupled so that the corresponding stress tensors are not separately conserved. In this setting, the ordinary charge conservation identities recover the factorization into two independent CFTs, whereas the pseudo-charge conservation identities hold for arbitrary values of the symmetry-breaking parameter. Finally, for conformal gravity in de Sitter space, we constrain the interactions between the graviton and the partially massless spin-2 field, both in the minimal theory and in the presence of additional scalar or vector matter, finding agreement with the predictions of the corresponding bulk theories.

We study the behaviour of general four-dimensional N=2 supergravity theories near boundaries at infinity of the vector multiplet moduli space, with a focus on a clear distinction between properties that follow from supergravity alone rather than from a consistent UV completion in quantum gravity. The basis for our analysis is the appearance of a real variation of polarised Hodge structure governing the vector multiplet moduli space of N=2 supergravity, combined with very recent advances in the mathematics literature. The resulting classification of infinite distance degenerations exactly agrees, at the level of gauge couplings, with a possible interpretation of the limits as decompactification limits to five or six dimensions or as weakly coupled string limits, as expected from a quantum gravity perspective. This observation extends previous results found in the context of string compactifications on Calabi–Yau varieties and, remarkably, holds for general four-dimensional N=2 supergravity theories regardless of a UV completion. We furthermore show the existence of a Kähler basis in which the triple couplings of the prepotential are non-negative, irrespective of a geometric origin of the supergravity theory. Finally we propose two non-trivial constraints for a four-dimensional N=2 supergravity to be compatible with quantum gravity: Integrality of the Hodge structure as a necessary condition for the completeness hypothesis and the existence of a simple normal crossing compactification of the moduli space to ensure compatibility of the weakly coupled asymptotic gravitational duality frames.

We explore a vast landscape of scale-covariant string backgrounds that arise from the backreaction of D$p$-branes. In a convenient Weyl frame the metric has an AdS factor and the dilaton has a non-trivial profile. This structure provides a notion of scale covariance that facilitates a holographic analysis and provides a window into the dynamics of the dual strongly coupled non-conformal gauge theory. We present different classes of supergravity solutions of this kind and discuss their dual field theory interpretation. We find several new supersymmetric backgrounds arising from D2-branes by utilizing four-dimensional gauged supergravity, D$p$-branes wrapped on compact manifolds, and infinite families constructed directly in ten dimensions. We construct holographic RG flows between some of the new supersymmetric backgrounds and compute the spectrum of low-lying local operators. Holography determines the form of observables such as the thermal free energy and scalar two-point functions, predictions which could be tested on the lattice.

Correlations in asymptotic fluxes provide one of the primary experimental means of studying physical systems. In the context of relativistic QFT, the formalization of such experiments in terms of detector operators has led to significant progress. Motivated by these successes, in this paper we initiate a study of detector operators in nonrelativistic conformal field theories (NRCFTs). In NRCFTs, fluxes are distributed in both angle and velocity, making the natural detector operator the energy flux differential in velocity, $\mathcal{E}_v(\hat{n})$. We provide non-perturbative definitions of these operators in NRCFTs, and study their symmetries and algebra. We calculate one-point functions of these detector observables in a variety of few and many body states of representative NRCFTs, including free fermions, as well as fermions at unitarity. In interacting three-body states, they exhibit a remarkably rich behavior reflective of the Efimov three-body wave-function, including the formation of an enhanced "jet" structure due to the presence of resonant interactions. We use these detectors to derive a nonrelativistic generalization of the Hofman-Maldacena conformal collider bounds by arguing for positivity of $\mathcal{E}_v(\hat{n})$, and imposing this in spin 1 states. We outline a number of future directions for the study of detectors in NRCFTs, and discuss how they can be measured in a variety of experimental platforms from cold atoms to low energy neutron experiments.

Recent progress in network renormalization has identified various ways to transform network representations consistently across resolution levels. In particular, the renormalization flow describing how random graph models transform under node aggregation has identified a fixed point corresponding to an aggregation-invariant model that preserves the form of the connection probability across arbitrary scales. While this probabilistic multiscale approach has successfully been used to model real-world networks, construct latent-space node embeddings, and reconstruct finer-grained versions of aggregate network data in a principled resolution-invariant way, it has so far been restricted to binary graphs. Here we reformulate multiscale network renormalization for weighted random graphs. Via a generating-function approach, we identify a weighted multiscale model whose full probability law, encoding both link creation and weight assignment, has an invariant compound-Poisson form under arbitrary node aggregation. The renormalization rules that transform parameters across aggregation levels allow the ensemble to be calibrated at one level and applied at another, enabling consistent modeling and reconstruction of weighted networks across arbitrary scales.

In this paper, an algorithm for constructing orbifolds of Calabi-Yau manifolds of the Berglund-Hubsch-Krawitz (BHK) type and their mirror pairs for loop-type polynomials is proposed and applied. A complete enumeration of loop polynomials satisfying the Calabi-Yau condition was carried out, yielding 216 unique polynomials. For each polynomial, the Euler characteristics, the numbers of generations and singlets were computed, and the twisted sector states were analyzed using Roan's generalized combinatorial construction. It was established that Roan pairs are absent for all studied configurations, indicating a zero contribution of the twisted sectors to the difference in Hodge numbers upon resolving the orbifold singularities. The primary physical result of this work is the discovery of two configurations with diagonal exponents $\{74, 3, 3, 5, 3\}$ and $\{74, 3, 5, 3, 3\}$, for which the number of quark and lepton generations is precisely three ($N_{\text{gen}} = 3$), which is of key interest for constructing phenomenologically realistic models of particle physics.

Using the effective field theory of inflationary fluctuations extended by a massive isocurvature scalar, at unit sound speed, we compute the primordial bispectra generated by every cubic interaction of dimension four or less in the unitary gauge, with the curvature–isocurvature quadratic mixing resummed to all orders in the dimensionless mixing strength $\lambda=\rho/H$. Placing the mixing in the free Hamiltonian [...], every scale-invariant tree-level bispectrum reduces to a single Schwinger-parameter integral over independent leg kernels, and what perturbation theory calls single, double and triple exchange are contact diagrams of the resummed theory. At strong mixing we unveil genuinely new, exact multifield bispectrum shapes with no correlation to the standard shapes used in data analyses. The squeezed limit carries a cosmological collider signal with a universal frequency set by the dressed mass parameter $\mu_{\rm eff}$ at every exchange order. At weak mixing the clock coefficients are obtained in closed form and reduce [...] to the known perturbative single- and double-exchange results, exactly; the triple-exchange leading coefficient is obtained in closed form at every mass, in terms of hypergeometric functions, for the first time. At strong mixing the collider amplitude scales as $e^{\pi(\lambda-2\mu_{\rm eff})/2}$ over the whole $(\lambda,\mu_{\rm eff})$ plane. This encompasses all previously studied regimes as particular cases [...]. A first confrontation of the six exact shapes with the Planck PR4 binned bispectrum shows a preference for strong mixing, where the cosmological collider oscillations invade mildly squeezed configurations [...]. The shapes of all interactions over the whole plane [...] are released with the code that produced them, so that a template at any mixing and mass can be read off in milliseconds.

Light scalar fields with gravitational strength interactions (moduli) are ubiquitous in string, supergravity, and extra-dimensional models. If produced in the early universe with lifetimes exceeding the epoch of recombination, slowly decaying moduli inject energetic particles into the Standard Model plasma, leaving distinct imprints on the cosmic microwave background (CMB). In this paper we systematically derive CMB injection bounds on light moduli ($10^3\ {\rm eV}\leq m_{\rm{mod}}\leq 10^8\ {\rm eV}$) across the primary cosmological production channels: vacuum misalignment, heavy scalar/inflaton decay, thermal emission from the primordial plasma including production from a Hagedorn phase of hot strings. We show that these injection bounds impose stringent constraints that are complementary to both the standard cosmological moduli problem and fifth force bounds.

We consider conformal boundaries and defects in Euclidean conformal field theory (CFT) at finite temperature and infinite volume. We distinguish defects that wrap the thermal circle and defects that are localised at a point on it. In both set-ups, and for any dimension and co-dimension, we classify bulk one-point functions of scalar, vector and spin-2 primaries. We relate the low-temperature expansion of bulk scalar one-point functions to zero-temperature defect CFT data and new finite-temperature data of defect primaries. For localised defects, we derive one-point bootstrap sum rules using periodicity and the defect operator expansion. We consider several examples in free bulk CFTs: unitary and non-unitary boundary CFTs, a wrapping monodromy defect, and a localised magnetic line. We determine thermal defect CFT data in all examples. For defects wrapping the thermal circle, we compute thermodynamic observables. For localised defects, we observe a Hurwitz zeta function decomposition of one-point functions suggestive of an underlying dispersion relation.

An emergent gauge field can control boundary critical behavior without changing the bulk quantum spin liquid. We investigate the boundary criticality in $(2+1)$-dimensional $U(1)$ Dirac spin liquid, whose low-energy physics is described by massless quantum electrodynamics. Using a perturbative renormalization-group analysis in a half-space, we show that Neumann and Dirichlet boundary conditions for the emergent gauge field lead to distinct boundary universality classes. We determine the boundary scaling dimensions of the fundamental fields and gauge-invariant operators that provide observable signatures. We further propose a fermion-gauge lattice model with tunable boundary interactions as a microscopic setting for realizing and probing these boundary universality classes.

We study the effective potentials for various probe D-branes surrounding AdS$_5$ black holes in consistent truncations of type IIB supergravity on Sasaki-Einstein 5-manifolds $Y^{p,q}$, focusing in particular on the conifold base $T^{1,1}$. The probes are either dual Giant Gravitons having support on AdS spacetime dimensions or D3 probe branes wrapped on non-contractible cycles of the internal manifold, which then appear "point-like" in the AdS spacetime. The non-contractible cycles of $T^{1,1}$ and $Y^{p,q}$ topologically stabilize the branes at a finite size, and we find stable supersymmetric wrapped D3 branes on Gutowski-Reall black holes, effectively providing a stable BPS composite system in AdS$_5$ made by a black hole surrounded by a baryonic condensate. We show that the onset of superradiance in the black hole backgrounds matches the emergence of stable minima in the probe effective potentials, extending the constructions of Dual Dressed Black Holes and the Grey Galaxies. We comment on the relevance of this result for recent studies concerning the classical cohomology problem and the associated fortuitous states in the dual Klebanov-Witten theory.

We study the inflationary no-boundary wavefunction in a “linear-roll limit,” in which $G_N \to 0$ while the Hubble scale $H$ and the slope of the inflaton potential are held fixed. Although bulk graviton fluctuations are suppressed in this limit, quantum inflaton fluctuations remain finite. Because reheating occurs on a surface of fixed inflaton value, these fluctuations become fluctuations of the conformal factor of the reheating metric. The probability measure for these fluctuations simplifies drastically and becomes Gaussian, even for large fluctuations. The resulting log-correlated theory at superhorizon scales is a $d$-dimensional analog of Liouville gravity without the exponential potential, which we call Logarithmic Gravity. The quadratic and linear terms are determined by conformally covariant scattering data; in even dimensions we relate them to the critical GJMS operator and Branson $Q$-curvature. From bulk unitarity, we argue that one-loop effects do not generate additional conformal-factor dependence in the physical probability measure, and we verify the cancellation explicitly in $2 + 1$ dimensions. For the no-boundary norm with a rolling inflaton, we show that the noncompact residual conformal symmetries can be gauge fixed with a finite Faddeev-Popov determinant, so they do not force the sphere contribution to vanish. In addition, conditioning on the absence of vacuum decay by bubble nucleation generates the Liouville exponential potential in the probability measure.

We investigate the recently introduced polygonal hypergeometric functions [arXiv:2502.12127, arXiv:2507.01904], which are conjectured to evaluate multipoint parametric one-loop conformal integrals in arbitrary dimensions. These functions can be systematically constructed using a diagrammatic algorithm that operates in terms of simple planar shapes. Given an $n$-point polygonal power series, we formulate the corresponding GKZ hypergeometric system, which admits an explicit diagrammatic interpretation. In particular, the transianic matrix, which encodes all the relevant diagrammatic data, determines the corresponding toric matrix along with a basis of the lattice of relations, thereby singling out a finite subsystem of $n(n-3)/2$ toric equations. This finite subsystem consists of second- and third-order equations, whereas the full toric family also includes higher-order ones. We illustrate the general procedure with several prominent examples of polygonal functions, including the fourth Appell and the Srivastava–Daoust functions.

We consider a planar and a non-planar two-loop Feynman integral family arising in non-factorisable quantum chromodynamics corrections to Higgs-boson production via vector boson fusion. We construct basis integrals which fulfil canonical differential equations with respect to seven variables for the external kinematics and an internal mass. One of the integral families involves nested square roots. We present strategies for deriving canonical master integrals, $\mathrm{d}\log$ forms and analytic continuation for such cases. Including also families with crossed massless legs, we express the $\varepsilon$-expanded integrals in terms of an algebraically independent function basis. For their evaluation, we employ numerical integration of the differential equations.

We formulate a class of systems of Euclidean correlation functions for compact $d$-dimensional conformal field theories which are not necessarily unitary. We also introduce conformal $d$-vertex algebras, $d$-dimensional analogues of full vertex algebras, and prove that such an algebra is obtained from every system of conformally covariant correlation functions. Finally, for each $d\geq 2$, we construct a family of examples satisfying these axioms, corresponding to generalized free scalar fields. For $d>2$, this family contains both unitary and non-unitary examples.

We study heavy-heavy-light-light correlators in holographic CFTs when the heavy operator is dual to a compact object localized on scales much smaller than the AdS radius. In the long-wavelength regime, the internal structure of the heavy state is described by a worldline effective field theory, whose local operators encode finite-size response coefficients such as tidal Love numbers. Focusing on weakly self-gravitating objects, we show that the Witten diagrams generated by local heavy-light interactions reorganize into the Born series of an effective wave equation for the light bulk field. We reproduce the same correlator by solving this wave equation with localized interactions and identify how the worldline Wilson coefficients are encoded in the CFT data. Our results provide a direct relation between the finite-size response of compact objects in AdS and the dimensions and OPE coefficients of heavy-light and light-light composite operators in the dual CFT.

We extend a lattice-QCD determination of the response of the quark-gluon plasma to an external magnetic field using a Taylor expansion of the pressure at zero field. A zero-net-flux “half-and-half” magnetic-field configuration permits continuous differentiation with respect to the field amplitude on gauge ensembles generated at $B=0$. The electromagnetic vector potential is chosen to minimize stochastic noise in the Taylor operators. Previous results through $O(B^2)$ are extended to $O(B^4)$ using $2+1$-flavor HISQ/tree ensembles with $m_l/m_s=0.05$ and $N_t=8$. The renormalized quadratic coefficient confirms paramagnetic behavior above the QCD crossover. The fourth-order coefficient $C_4$ is small at high temperature and is enhanced near the crossover. At $T=167$ MeV the $C_4$ result at $48^3\times 8$ lies about two standard deviations above the corresponding $32^3\times8$ one, indicating likely finite- volume effects. For $eB = 0.2$ GeV$^2$, characteristic of magnetic fields relevant to heavy-ion phenomenology, the fourth-order contribution is small compared with the leading magnetic correction at moderate and high temperatures, but near the crossover the large-volume result is about half of the central value of the quadratic term, after the vacuum-pressure subtraction. We also place the method in the context of other lattice determinations of magnetic susceptibility and clarify the relevance of finite-volume effects for the thermal subtraction. This work was completed in 2014 and preliminary results were reported at Lattice 2014.