# A relativistic continuous matrix product state study of field theories with defects

Karan Tiwana, Edoardo Lauria, Antoine Tilloy

► To cite this version:

Karan Tiwana, Edoardo Lauria, Antoine Tilloy. A relativistic continuous matrix product state study of field theories with defects. Journal of High Energy Physics, 2025, 2025 (05), pp.097. 10.1007/JHEP05(2025)097. hal-04919697

HAL Id: hal-04919697

https://hal.science/hal-04919697v1

Submitted on 27 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# A relativistic continuous matrix product state study of field theories with defects

Karan Tiwana<sup>1, ∗</sup>, Edoardo Lauria<sup>1</sup> and Antoine Tilloy<sup>1</sup>

Laboratoire de Physique de l'École Normale Supérieure, Mines Paris - PSL, Inria, CNRS,
ENS-PSL, Sorbonne Université, PSL Research University,
Paris, France

E-mail: karantiwana21693@gmail.com, edoardo.lauria@minesparis.psl.eu,
antoine.tilloy@minesparis.psl.eu

**ABSTRACT:** We propose a method to compute expectation values in 1+1-dimensional massive Quantum Field Theories (QFTs) with line defects using Relativistic Continuous Matrix Product State (RCMPS). Exploiting Euclidean invariance, we use a quantization scheme where (imaginary) time runs perpendicularly to the defect. With this choice, correlation functions of local operators in the presence of the defect can be computed as expectation values of extended operators in the no-defect vacuum, which can be approximated by a homogeneous RCMPS. We demonstrate the effectiveness of this machinery by computing correlation functions of local *bulk* and *defect* operators in $\phi^4$ theory with a magnetic line defect, in perturbative, strong coupling, critical, and symmetry-broken regimes.

**KEYWORDS:** Boundary Quantum Field Theory, Field Theories in Lower Dimensions

ARXIV ePRINT: 2501.09797
# Contents

<table><tbody><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>1.1</td><td>Defects</td><td>2</td></tr><tr><td>1.2</td><td>The variational method</td><td>3</td></tr><tr><td>1.3</td><td>Rotating the defect, from impurity to extended operator</td><td>3</td></tr><tr><td>1.4</td><td>Advantages and limitations</td><td>4</td></tr><tr><td>2</td><td>Setup</td><td>5</td></tr><tr><td>2.1</td><td>A scalar field model without defects</td><td>5</td></tr><tr><td>2.2</td><td>Defects</td><td>6</td></tr><tr><td>3</td><td>Relativistic continuous matrix product states</td><td>7</td></tr><tr><td>3.1</td><td>A motivation</td><td>7</td></tr><tr><td>3.2</td><td>The definition of RCMPS</td><td>8</td></tr><tr><td>3.3</td><td>Computing local observables</td><td>9</td></tr><tr><td>3.4</td><td>Optimizing the state</td><td>11</td></tr><tr><td>3.5</td><td>Introducing defects</td><td>11</td></tr><tr><td>4</td><td>Results</td><td>13</td></tr><tr><td>4.1</td><td>Perturbation theory calculations</td><td>13</td></tr><tr><td>4.2</td><td>Weak coupling g = 0.1, comparison with perturbation theory</td><td>15</td></tr><tr><td>4.3</td><td>Strong coupling g = 2.0 self-comparison</td><td>16</td></tr><tr><td>4.4</td><td>Critical coupling g ≈ 2.771525, comparison with scaling expectations</td><td>17</td></tr><tr><td>4.5</td><td>Symmetry broken phase, coupling g = 4</td><td>19</td></tr><tr><td>5</td><td>Conclusion</td><td>21</td></tr><tr><td>A</td><td>Derivation of RCMPS formulae for defects</td><td>22</td></tr><tr><td>A.1</td><td>Expectation value of the defect operator</td><td>23</td></tr><tr><td>A.2</td><td>One-point functions in the defect model</td><td>24</td></tr><tr><td>A.3</td><td>Strategy for more general expectation values</td><td>26</td></tr><tr><td>B</td><td>No-defect two-points correlators</td><td>27</td></tr><tr><td>B.1</td><td>Perturbative regime</td><td>27</td></tr><tr><td>B.2</td><td>Critical coupling regime</td><td>28</td></tr></tbody></table>

## 1 Introduction

Relativistic scalar quantum field theories (QFT) in 1 + 1 space-time dimensions are an excellent sandbox for theorists: they can be defined rigorously in a fairly straightforward way, and yet are generically difficult to solve non-perturbatively. Intuitively, such theories are numerically challenging because their number of degrees of freedom grows both when the
short distance (UV) cutoff and long distance (IR) cutoff are lifted. Nonetheless, a tremendous amount of progress has been made in recent years to compute local expectation values of low dimensional QFT beyond perturbation theory. In this article, our objective is to go beyond simple local expectation values, and get numerically precise predictions in 1+1 dimensional QFTs with defects. This has mainly two interests: 1) we get numerical access to new non-trivial physics, like that of impurities, and 2) we can benchmark non-perturbative methods on a more difficult example.

## 1.1 Defects

In the general case, i.e. in a $d$-dimensional space-time, a $p$-dimensional defect is a modification of the $d$-dimensional action (the “bulk theory”) that has support on $p$-dimensional submanifold $B$ of $\mathbb{R}^d$. More explicitly, but at a (so far) non-rigorous level, we consider “bulk” QFTs with an action of the form

$$S_{\text{bulk}}(\phi) = \int_{\mathbb{R}^d} \frac{(\nabla \phi)^2}{2} + V(\phi), \qquad (1.1)$$

where we will later fix $d = 2$. The “defect” QFT is specified by the action

$$S(\phi) = S_{\text{bulk}}(\phi) + S_{\text{defect}}(\phi) = S_{\text{bulk}}(\phi) + \int_B f(\phi). \qquad (1.2)$$

Expectation values can be computed with this modified action in the same way as for a standard QFT. Again non-rigorously, in Euclidean signature we have that:

$$\langle \phi(x_1) \cdots \phi(x_n) \rangle_{\text{defect}} = \frac{\int \mathcal{D}\phi \, \phi(x_1) \cdots \phi(x_n) e^{-S(\phi)}}{\int \mathcal{D}\phi \, e^{-S(\phi)}} = \frac{\langle \phi(x_1) \cdots \phi(x_n) e^{-\int_B f(\phi)} \rangle_{\text{bulk}}}{\langle e^{-\int_B f(\phi)} \rangle_{\text{bulk}}}. \qquad (1.3)$$

Physically, defects can be used to model boundaries, interfaces, and impurities which are very common in realistic low-energy systems (e.g. in the Kondo effect [1, 2] or surface criticality, see [3, 4] for review). Considering defects is also useful in the formal study of QFT, and e.g. allows for a characterization of symmetries that are not captured by local QFT operators [5] (see [6] for a recent review).

There has been a lot of progress in understanding the role of (non-topological) defects in QFT, mostly when the bulk action QFT is conformally invariant, both in $d = 2$ [7–9] and higher (see e.g. [10, 11] for early studies) with a recent explosion [12–27]. However, much less is known for gapped bulk actions (conformally invariant at short but not large distance) beyond perturbative regimes.

In this work, we study a class of (non-topological) line defects for the massive $\phi^4$ theory: the so called magnetic line defects. As reviewed in [28–33] and discussed more recently in [13, 20], such defects play an important role in condensed matter systems where they capture the physics of magnetic-like impurities. For conformal QFTs in $d \ge 2$, magnetic line defects can be studied with a variety of techniques, starting from perturbation theory and conformal bootstrap [14–18, 34–41] and more recently fuzzy sphere regularization [27].
## 1.2 The variational method

Independently of defects, the two historical methods to solve generic QFTs are perturbation theory and lattice Monte Carlo, both of which have limitations. In the first case, one can lift both UV and IR cutoffs, but results are given as a divergent power series in the coupling (which can be Borel resummed in favorable cases, see e.g. [42, 43]). In the second case, results are fully non-perturbative, at the cost of hard UV and IR cutoffs (a finite size space-time lattice), and statistical errors. There is now a plethora of less traditional methods, including the modern S-matrix bootstrap (see [44] for a recent review and references), functional renormalization [45, 46], and possibly soon even direct quantum simulation [47–49], that could in principle be modified to account for defects.

The variational approach yields another class of non-perturbative methods (including Hamiltonian truncation [50–53]) which we now focus on. The idea is fairly simple: treat a QFT in canonical quantization simply as a particular quantum mechanical system, with a Hilbert space $\mathcal{H}$ and Hamiltonian

$$H_{\text{bulk}} = \int_{\mathbb{R}} h := \int_{\mathbb{R}} \frac{\hat{\pi}^2}{2} + \frac{(\nabla \hat{\phi})^2}{2} + V(\hat{\phi}), \quad \text{with} \quad [\hat{\phi}(x), \hat{\pi}(y)] = i\delta(x-y), \quad (1.4)$$

and approximate its ground state (the non-perturbative vacuum) $|\Omega\rangle$ by minimizing the energy density over a manageable submanifold $\mathcal{M}$ of states:

$$|\Omega\rangle \approx |\psi_0\rangle := \underset{|\psi\rangle \in \mathcal{M} \subset \mathcal{H}}{\operatorname{argmin}} \frac{\langle\psi|h\rangle}{\langle\psi|\psi\rangle}. \quad (1.5)$$

The output of a bona fide variational approach is an approximate ground state $|\psi_0\rangle$ and an efficient algorithm to evaluate expectation values $\langle\psi_0|\hat{\phi}(x_1)\cdots\hat{\phi}(x_n)|\psi_0\rangle$ of interest. This strategy is particularly adapted to 1 + 1 dimensional scalar QFT because most models have a finite energy density in the vacuum, and thus the optimization problem is in principle well posed.

A particularly powerful class of states to use in the variational approach are relativistic continuous matrix product states (RCMPS) [54, 55]. They have a number of appealing features: 1) they are defined directly in the continuum (without UV cutoff), 2) their thermodynamic limit can be taken easily (no IR cutoff), 3) their expressiveness (the dimension of $\mathcal{M}$) can be arbitrarily increased, 4) local expectation values can be computed (numerically) exactly, and 5) optimizing over them is practically efficient. These states have been successfully used to study vacuum expectation values in the $\phi^4$ [55], Sine-Gordon, and Sinh-Gordon models [56].

## 1.3 Rotating the defect, from impurity to extended operator

In $d = 1 + 1$ space-time dimension, we can consider line defects (i.e. $p = 1$): a one dimensional deformation of the two dimensional action, which we take to be a straight line. In the Hamiltonian picture, such a defect is naturally interpreted as an impurity, namely a perturbation of the Hamiltonian localized at one point $x_0$ in space, and persistent in time:

$$H = H_{\text{bulk}} + H_{\text{defect}} = H_{\text{bulk}} + f[\hat{\phi}(x_0)]. \quad (1.6)$$

If the defect is only a finite or semi-infinite line, this corresponds to a local quench, that is an impurity that is turned on only for some time.
<figure><img src="image_3.png" /><figcaption>Figure 1. a) original defect, interpreted as an impurity, b) equivalent rotated defect, interpreted as an extended operator c) representation of the defect expectation value in the Hamiltonian formalism: the bulk dynamics from $+\infty$ and $-\infty$ is equivalent to a projection of the state at fixed $\tau$ to the bulk ground state $|\Omega\rangle$.</figcaption></figure>

For generic quantum many-body systems, such impurities/quenches give rise to rich physics but introduce challenges with the variational method. First, translation invariance in space is lost. Second, in the local quench case, the quantum state needed to evaluate expectation values becomes time dependent, and thus needs to be evolved e.g. with the time dependent variational principle (TDVP) [57, 58].

Fortunately, a simple observation makes a plain variational study possible for large subset of defects and operator insertions. For a relativistic QFT in imaginary time, i.e. a Euclidean invariant theory, we can rotate the quantization direction, exchange the role of space and time, and make the defect spacelike instead of timelike (see figure 1). In this new picture, time evolution is generated by the bulk Hamiltonian, and expectation values are simply evaluated by applying an extended operator to the translation-invariant bulk vacuum, $|\Omega\rangle$:

$$ \langle\phi(x_1) \cdots \phi(x_n)\rangle_{\text{defect}} = \frac{\langle\Omega| \hat{\phi}(x_1) \cdots \hat{\phi}(x_n) e^{-\int_B f(\hat{\phi})} |\Omega\rangle}{\langle\Omega| e^{-\int_B f(\hat{\phi})} |\Omega\rangle} \quad (1.7) $$

This representation is convenient because $|\Omega\rangle$ can be replaced by its variational approximation as a RCMPS $|\psi_0\rangle$, simply computed from the bulk theory. The last ingredient is technical, non-trivial, and specific to RCMPS: if the operator insertions and defect are aligned (same imaginary time), and if $f$ is linear (i.e. a so called magnetic defect), the right-hand side of (1.7) is efficiently computable at the same asymptotic cost with the defect than in the vanilla no defect case. All the complex optimization toolbox required to find RCMPS ground states for translation invariant models can thus be reused without modification. All that is needed is a particular post-processing of this ground state, which is comparatively cheap.

## 1.4 Advantages and limitations

Although we test it only for a $\phi^4$ potential here, the method can be used for any bulk model with a potential $V(\phi)$ for which RCMPS optimization algorithms have been developed. Currently, this includes all polynomial potentials, real exponentials $e^{\beta\phi}$, and imaginary exponentials $e^{i\beta\phi}$ with $\beta < \sqrt{4\pi}$.

The method we present has the same cost and precision as standard RCMPS. At short distances, in the deep UV, RCMPS behave like the free Fock vacuum. For the super-renormalizable bulk theories we consider, this is also the behavior of the true ground state,
<figure><img src="image_4.png" /><figcaption>**Figure 2.** a) an allowed configuration of defect (continuous line) and local operators (cross) that we can compute because all operator insertions are aligned b) a configuration that a plain variational method cannot deal with, because the operator insertions are not aligned with the defect. The latter setup could be accessed, in principle, with an evolution method like TDVP.</figcaption></figure>

and thus our method has excellent precision in this limit. As we will see, it compares favorably with third-order perturbation theory. At large distance, RCMPS provide good approximations of gapped theories even at strong coupling, and we will observe that this precision is transferred to the defect case. Finally, if the bulk is fine-tuned to criticality, RCMPS no longer provide as good an approximation of long distance physics, and introduce an effective length-scale associated to the limited amount of entanglement they can represent. Precision is then reduced and extrapolations to the infinite entanglement limit [59, 60] would be required. This is the first limitation of our method which makes it complementary to previous studies that considered scale invariant field theories in the bulk.

While the bulk model can be very general, we are limited in the type of defects we can study efficiently. First, the defects have to be linear in the field $\phi$, i.e. correspond to operator insertions of the form $\exp(-\mu \int_0^L \phi)$. The geometry has flexibility: we allow combinations of defects of arbitrary length $L$ (even numerically semi-infinite or infinite) and strength $\mu$, along with arbitrary local operator insertions, as long as they are all aligned (see figure 2 for allowed configurations). This latter point is the second limitation of the method, and likely of any vanilla variational approach. This nonetheless includes impurity and complex local quench problems, as long as one is interested only in measurements on the impurity. In the local quench case, we remain limited to imaginary time, which is the third limitation of our approach at its present level of development.

## 2 Setup

In this section, we explain how the $\phi^4$ model can be defined easily in the operator representation and then introduce a magnetic defect.

### 2.1 A scalar field model without defects

The $\phi^4$ model is (so far naively) defined by the Hamiltonian $H_{\text{bulk}}$

$$H_{\text{bulk}} := \int_{\mathbb{R}} \frac{\hat{\pi}^2}{2} + \frac{(\nabla \hat{\phi})^2}{2} + \frac{m^2}{2} \hat{\phi}^2 + g \hat{\phi}^4 = H_m + \int_{\mathbb{R}} g \hat{\phi}^4, \quad (2.1)$$

where $m$ is the bare mass and $g$ is the coupling constant. At this stage, this Hamiltonian is not yet well defined and one option to make sense of it is to first diagonalize the free
part $H_m$ in terms normal modes $\hat{a}_k$. More explicitly, we expand the field $\hat{\phi}(x)$ and its conjugate momentum $\hat{\pi}(x)$

$$ \hat{\phi}(x) = \frac{1}{2\pi} \int \frac{dk}{\sqrt{2\omega_k}} (\hat{a}_k e^{ikx} + \hat{a}_k^\dagger e^{-ikx}) \qquad (2.2) $$

$$ \hat{\pi}(x) = \frac{1}{2i\pi} \int dk \sqrt{\frac{\omega_k}{2}} (\hat{a}_k e^{ikx} - \hat{a}_k^\dagger e^{-ikx}), \quad \text{with} \quad [\hat{a}_k, \hat{a}_{k'}] = 2\pi\delta(k-k'), \qquad (2.3) $$

where $\hat{a}_k^\dagger$ and $\hat{a}_k$ are the creation and annihilation operators (respectively) of mode $k$, and $\omega_k = \sqrt{k^2 + m^2}$. This yields:

$$ H_m = \int_{\mathbb{R}} \frac{dk}{2\pi} \frac{\omega_k}{2} (\hat{a}_k \hat{a}_k^\dagger + \hat{a}_k^\dagger \hat{a}_k). \qquad (2.4) $$

The Hamiltonian above describes infinitely-many decoupled harmonic oscillators, one for every value of $k$. The Hilbert space of this model is constructed by adding the excitations on top of the free vacuum $|0\rangle$ which is annihilated by $\hat{a}_k$'s. Normal ordering the free Hamiltonian $H_m$ w.r.t. mass $m$ sets the ground state energy density to zero. Now we have a legitimate Hamiltonian which is bounded from below, albeit without any interactions:

$$ :H_m: := \int_{\mathbb{R}} \frac{dk}{2\pi} \omega_k \hat{a}_k^\dagger \hat{a}_k. \qquad (2.5) $$

We can further apply this normal ordering to the interaction potential. In terms of Feynman diagrams, this normal ordering is equivalent to the suppression of tadpole graphs, which are the only divergences that appear in perturbation theory. In $d=1+1$, it has been shown in [61, 62] that normal ordering polynomial interactions is indeed sufficient beyond perturbation theory, and leads to a Hamiltonian density that is bounded from below (at least as long as monomial of largest degree is even and with a positive coupling). Therefore, the $\phi^4$ Hamiltonian

$$ H_{\text{bulk}} = :H_m: + \int_{\mathbb{R}} g :\hat{\phi}^4: = \int_{\mathbb{R}} h, \qquad (2.6) $$

is amenable to the variational method: in the translation invariant case, we minimize a finite, lower-bounded scalar to get $\langle \Omega|h|\Omega \rangle$.

This model is arguably the simplest interacting QFT that is rigorously defined, while not being supersymmetric or integrable. As a result, it has been studied with various methods, each with its strengths and weaknesses: Hamiltonian truncation [51–53] (IR and energy cutoff), Monte-Carlo [63–65] (both UV and IR cutoffs), MPS [66, 67] (UV cutoff), Tensor Network Renormalization [68, 69] (UV cutoff), Borel Resummation [42] (no cutoffs but perturbative). Aside from not being limited by any cutoff, RCMPS can straightforwardly handle any interaction potential added to $H_m$ as long as the Hamiltonian density is bounded from below, with the same construction that is used for the $\hat{\phi}^4$ model. In this work, we restrict ourselves to the $\hat{\phi}^4$ model for simplicity and concreteness.

## 2.2 Defects

Next, we introduce a magnetic line defect of length $L$ in the model. In the Hamiltonian picture, this is a modification of the $\hat{\phi}^4$ Hamiltonian (2.6) where we turn on the additional
interaction $-\mu\hat{\phi}(x_0)$ localized at $x_0 = (x, 0)$, for an interval of imaginary-time $x$ of length $L$. As explained in 1.3, upon exploiting Euclidean invariance of the $\phi^4$ model, we can view such modification as an extended operator

$$ \mathcal{D}_L := e^{-\mu \int_{-L}^{0} \hat{\phi}}, \qquad (2.7) $$

so that vacuum expectation values in this model can be computed as

$$ \langle\phi(x_1) \cdots \phi(x_n)\rangle_{\text{defect}} = \frac{\langle 0, g | \hat{\phi}(x_1) \cdots \hat{\phi}(x_n) \mathcal{D}_L | 0, g \rangle}{\langle 0, g | \mathcal{D}_L | 0, g \rangle}, \qquad (2.8) $$

where $|0, g\rangle$ is the interacting vacuum of $H_{\text{bulk}}$ at coupling $g$. The defect expectation value is both UV and IR finite<fnref n="1" /> without normal ordering.

## 3 Relativistic continuous matrix product states

In this section we motivate the definition of RCMPS, explain how to compute the expectation value of local observables, quickly discuss how RCMPS are optimized to find the ground state of a quantum field theory, and finally show how extended operators (defects) can be evaluated. Apart from the last part about defects, this section is a summary of previous results [54–56].

### 3.1 A motivation

Matrix product states are ansatz quantum states for spin chains (or chains of qudits). They are parameterized non-linearly and in an extensive manner. Moreover, for the translation invariant case, the number of parameters does not grow with system size. They efficiently approximate states with a bounded amount of bipartite entanglement, such as the ground states of local gapped Hamiltonians. A natural strategy to solve a QFT like $\phi^4$ is thus to discretize it on a chain with lattice spacing $\epsilon$, solve it with MPS, and finally numerically take the continuum limit $\epsilon \to 0$ [66]. This works fairly well practically for fixed $\epsilon$, but one cannot take $\epsilon$ arbitrarily small. Indeed, a relativistic QFT, even massive, behaves like a critical theory at short distance. As a consequence, the amount of bipartite entanglement in the ground state of the lattice discretized model grows logarithmically with $\epsilon^{-1}$. This makes the MPS approximation worse and worse as the continuum limit is approached (for a fixed number of variational parameters). Taking the continuum limit of MPS first analytically, with CMPS, does not bypass the divergent entanglement entropy problem. So while MPS do solve the IR problem admirably, they struggle with the short distance behavior that is typically found in relativistic QFT.

Hamiltonian truncation, on the other hand, has almost the opposite qualities. The idea is to first introduce a box of size $L$ (or equivalently an infrared cutoff), which has the effect of discretizing the momenta $k$ appearing in the normal modes $\hat{a}_k$. The resulting Hilbert space, spanned by all vectors of the form:

$$ |k_1, m_1, k_2, m_2, \cdots k_n, m_n\rangle := \frac{1}{\sqrt{m_1! m_2! \cdots m_n!}} a_{k_1}^{\dagger m_1} a_{k_2}^{\dagger m_2} \cdots a_{k_n}^{\dagger m_n} |0\rangle, \qquad (3.1) $$

<footnote n="1">¹By IR finite, we mean that the expectation value has a well defined limit when the defect size is fixed, but the thermodynamic limit is taken for the bulk, that is without global IR cutoff.</footnote>
where $|0\rangle$ is the free vacuum, is still an infinite dimensional Fock space, but with a discrete basis. One can further truncate it by considering only basis states below a certain cutoff energy $E_c$. Ultimately, one diagonalizes the Hamiltonian on the resulting finite dimensional Hilbert space to obtain the ground state. As far as short distance physics is concerned, this is an excellent scheme because the Hilbert space is constructed from excitations above $|0\rangle$. The latter, for super-renormalizable models, already contains an infinite amount of spatial entanglement, that diverges in exactly the same way as that of the true interacting ground state. However, just like exact diagonalization on the lattice, Hamiltonian truncation does not solve the many-body problem (or IR problem): for a fixed target error, the size of the truncated Hilbert space increases exponentially with system size. This is because, in contrast with matrix product state, the parameterization on the truncated Fock space is not extensive.

Ideally, one would like to benefit from the progress made both with HT and MPS and solve IR and UV difficulties together. This is doable: one should just construct a (continuous) matrix product state on top of the free Fock vacuum. This combines the insight from HT that one should start from a spatially entangled state with the same (real space) divergence structure as the true ground state, with the intuition that MPS provide a good extensive parameterization adapted to the thermodynamic limit.

## 3.2 The definition of RCMPS

To construct a (continuous) MPS, one first needs to write the Hilbert space as a (continuous) tensor product of factors. Choosing factors associated to each normal mode would be natural, but not adapted because 1) physics is not translation invariant in momentum 2) the interacting ground state is not weakly entangled in momentum. Following [54, 55], we go back to real space and introduce $\hat{a}^\dagger(x)$ which is the Fourier transform of the canonical creation operators $\hat{a}_k^\dagger$:

$$ \hat{a}(x) = \frac{1}{2\pi} \int_{\mathbb{R}} dk \, e^{ikx} \hat{a}_k. \qquad (3.2) $$

Note the absence of $\sqrt{\omega_k}$ factors in the denominator which are crucial to ensure that $a(x)$ has canonical commutation relation $[a(x), a^\dagger(y)] = \delta(x - y)$. This makes the Hilbert space a continuously infinite tensor product of factors associated with each oscillator at position $x$. This is all that is needed to define a CMPS.

Formally, a RCMPS on a line of size $L$ and with periodic boundary conditions is a quantum state parameterized by two $D \times D$ complex matrices $Q, R$ and defined as

$$ |Q, R\rangle = \text{tr}\{\mathcal{P} \exp \left[ \int_{0}^{L} dx Q \otimes 1 + R \otimes \hat{a}^{\dagger}(x) \right] \}|0\rangle, \qquad (3.3) $$

where $|0\rangle$ is the free Fock vacuum verifying $\forall x, \hat{a}(x)|0\rangle = 0$, $\mathcal{P} \exp$ is the path-ordered exponential, and the trace is taken over the finite $D$-dimensional space on which $Q, R$ act. This representation (3.3) makes the extensiveness of the state manifest, but it helps to
expand<fnref n="2" /> it explicitly in the Fock basis:

$$|Q, R\rangle = \sum_{n=0}^{+\infty} \int_{0 \le x_1 \le \dots \le x_n \le L} dx_1 \dots dx_n \varphi_n(x_1, x_2, \dots, x_n) a^\dagger(x_1) \dots a^\dagger(x_n) |0\rangle \quad (3.4)$$

$$\text{with } \varphi_n(x_1, x_2, \dots, x_n) := \text{tr} \left[ e^{x_1 Q} R e^{(x_2-x_1)Q} R \dots e^{(x_n-x_{n-1})Q} R e^{(L-x_n)Q} \right]. \quad (3.5)$$

This shows that while the RCMPS has few parameters, it is *not* sparse in the Fock space, and all the basis coefficients are generically non-zero. The reader unfamiliar with MPS techniques may very well forget the motivation 3.1 of RCMPS and take its expanded form (3.4) as starting point.

We now make a few remarks about the RCMPS we just defined. First, the parameterization with $Q$ and $R$ is redundant. The state $|Q, R\rangle$ is invariant under $Q \to X^{-1}QX$ and $R \to X^{-1}RX$. We can exploit this freedom to impose a gauge condition $Q + Q^\dagger + R^\dagger R = 0$, the so-called *left canonical gauge*. With this choice, which we make from now on, one can write $Q = -iK - R^\dagger R/2$ with $K$ self-adjoint. Second, and with this choice of gauge, one can take the thermodynamic limit seamlessly, because the state has constant norm in the infinite size limit. Thus, at least within correlation functions, one can replace the integration interval $[0, L]$ in (3.3) by $\mathbb{R}$.

Finally, we note that the field operator $\hat{\phi}(x)$ is not a local function of $a(x)$ and $a^\dagger(x)$, but is rather obtained from a convolution (see (3.8) below). This is where RCMPS differ from standard CMPS, where one would take oscillator modes $\hat{\psi}(x)$ local in the field, i.e. with $\hat{\psi}^\dagger(x) \propto \sqrt{\Lambda}\hat{\phi}(x) - i\hat{\pi}(x)/\sqrt{\Lambda}$ with some arbitrary UV scale $\Lambda$, and replace $|0\rangle$ by the state $|0\rangle_\psi$ annihilated by all the $\psi(x)$. While the Hamiltonian is strictly local in this representation, its ground state is infinitely entangled. RCMPS use a tensor factorization that breaks the strict locality of $H$, but yields a finitely entangled ground state.

### 3.3 Computing local observables

In many-body systems, and even more so in field theories, knowing a state does not necessarily imply that one can compute observables efficiently. We now show that, fortunately, we do get local observables at a cost $\propto D^3$ for RCMPS.

To this end, first note that since $\hat{a}^\dagger(x)$ satisfies the same commutation relation as $\hat{\psi}^\dagger(x)$, the results for expectation values for strings of $\hat{a}(x)$ follow from “standard” CMPS computations.

In particular, normal-ordered correlation functions of products of $a$ and $a^\dagger$ can be computed by taking derivatives of the generating functional

$$Z_{j',j} := \frac{\langle Q, R | \exp[\int j'(x) \hat{a}^\dagger(x)] \exp[\int j(y) \hat{a}(y)] | Q, R \rangle}{\langle Q, R | Q, R \rangle}, \quad (3.6)$$

with respect to the sources $j(x)$ and $j'(y)$. The generating functional admits a compact form in terms of $T, j(x), j'(y)$ [70]:

$$Z_{j',j} = \frac{\text{tr}[\mathcal{P} \exp(\int T + jR \otimes 1 + j'1 \otimes R^*)]}{\text{tr}[\mathcal{P} \exp(\int T)]}, \quad (3.7)$$

<footnote n="2">²This expansion is done e.g. in [70]. One option to obtain the result in a simple way is to introduce $R(x) := e^{Qx} R e^{-Qx}$, write $|Q, R\rangle$ as $\text{tr}\{\mathcal{P} \exp[\int_0^L dx R(x) \hat{a}^\dagger(x)]\}|0\rangle$ and then Taylor expand the path-ordered exponential.</footnote>
where from now on we may assume the integrals are over $\mathbb{R}$. This is the central formula for CMPS from which every correlation function can in principle be derived.

Of course, we are interested in correlation functions of local operators such as $\hat{\phi}(x)$, not $\hat{a}(x)$. However, the latter can be obtained from those of $\hat{a}(x)$ and $\hat{a}^{\dagger}(x)$ through the convolution

$$ \hat{\phi}(x) = \int dy J(x-y)[\hat{a}(y) + \hat{a}^{\dagger}(y)], \quad J(x) = \frac{1}{2\pi} \int \frac{dk}{\sqrt{2\omega_k}} e^{-ikx}. \qquad (3.8) $$

This convolution is not too non-local: $J(x)$ has an integrable singularity at $x = 0$, but decays exponentially as $|x| \to \infty$. It only introduces a mild inconvenience compared to the CMPS scenario but is not inhibitory in practice due to the fast decay of $J(x)$. In fact, it is precisely because we already pay the price of some non-locality in the RCMPS case that computing extended observables (or equivalently operators in the complete defect model) comes with no extra cost as we discuss in 3.5.

In principle, one could compute expectation values for polynomials of $: \hat{\phi}(x) :$, by first computing expectation values of $\hat{a}(x)$ and then numerically convolving them with $J$. This would become increasingly costly as the degree of the polynomial increases. Fortunately, there is an easier way, which is to first compute the expectation value of a vertex operator $V_b(x) =: e^{b\hat{\phi}(x)} :$. The latter is indeed very simple in terms of the RCMPS generating functional $Z_{j',bJ}$,

$$ \langle V_b(0) \rangle = Z_{bJ,bJ}. \qquad (3.9) $$

Now, we may obtain expectation values of powers of the field by differentiating $\langle V_b(x) \rangle$ w.r.t. $b$ and taking $b = 0$ at the end. To compute $Z_{bJ,bJ}$, we simply interpret the path ordered exponential as the solution of an ordinary differential equation (ODE):

$$ \mathcal{O}(L) := \mathcal{P} \exp \left( \int_0^L \mathcal{T} + jR \otimes 1 + j'1 \otimes R^* \right) \quad \text{is indeed implied by} \qquad (3.10) $$

$$ \frac{d}{dx} \mathcal{O}(x) = \mathcal{O}(x) (\mathcal{T} + jR \otimes 1 + j'1 \otimes R^*) \quad \text{and} \quad \mathcal{O}(0) = 1. \qquad (3.11) $$

At this stage, it is customary to move to a superoperator representation $\mathcal{T} \mapsto \mathcal{L}$ to make explicit the fact that applying the map $\mathcal{T}$ scales as $D^3$ and not $D^4$ as it would naively if $\mathcal{T}$ were a generic map. Concretely, the superoperator representation is defined via the isomorphism $|w\rangle \otimes |v\rangle \cong |w\rangle\langle v|$ i.e. $\mathcal{L}$ is the representation of $\mathcal{T}$ acting on $D \times D$ matrices instead of vectors of lengths $D^2$. Explicitly:

$$ \mathcal{T}(|w\rangle \otimes |v\rangle) = Q|w\rangle \otimes 1|v\rangle + 1|w\rangle \otimes Q^*|v\rangle + R|w\rangle \otimes R^*|v\rangle \qquad (3.12) $$

$$ \Rightarrow \mathcal{L}(|w\rangle\langle v|) = Q(|w\rangle\langle v|) + (|w\rangle\langle v|)Q^{\dagger} + R(|w\rangle\langle v|)R^{\dagger}. \qquad (3.13) $$

Hence, for a generic matrix $\rho$,

$$ \mathcal{L} \cdot \rho = Q \rho + \rho Q^{\dagger} + R \rho R^{\dagger}. \qquad (3.14) $$

As a result, applying $\mathcal{L}$ can be done with simple matrix multiplications at cost $D^3$. Note that with the left-canonical gauge choice we made before, $\mathcal{L}$ is of the Lindblad form (and is thus in particular trace preserving).
Now combining the ODE and superoperator representations, we have $\langle V_b \rangle = \text{tr}(\rho_\infty)$ with $\rho_\infty = \lim_{x \to \infty} \rho(x)$,

$$ \frac{d}{dx} \rho(x) = \mathcal{L} \cdot \rho(x) + bJ(x)[R\rho(x) + \rho(x)R^\dagger], \qquad (3.15) $$

and initial condition $\rho(-\infty) = \rho_0$, with $\mathcal{L} \cdot \rho_0 = 0$ and $\text{tr}[\rho_0] = 1$.

To compute expectation values of powers of the field, we may now *forward differentiate* this ordinary differential equation (3.15) and evaluate the result at $b = 0$. The kinetic part of the Hamiltonian density (2.1) is slightly more complicated to evaluate, since it also involves derivatives of the field and the conjugate momentum, but it can be computed with similar techniques and the same asymptotic cost [54].

## 3.4 Optimizing the state

We may therefore obtain the expectation value of the Hamiltonian density, $\langle h \rangle$, as the trace of a solution of linear matrix ODEs involving $Q$ and $R$. These ODEs can be evaluated to essentially arbitrary numerical precision using high-order Runge-Kutta solvers, with a cost dominated by the application of the generator, hence $D^3$. Further, using backpropagation (or equivalently adjoint methods) we can compute the full gradient of $\langle h \rangle$ with respect to the relevant parameters $K, R$ at the same cost $D^3$.

We have everything needed to minimize the energy density over the parameters, and thus obtain a variational approximation to the ground state. The minimization could a priori be done using any numerical solver (plain gradient descent or LBFGS). However, as expected from conventional tensor network wisdom, it was observed [54–56] that such naive minimizations work poorly, and get stuck in plateaus. The solution is to incorporate the geometry of RCMPS, that is, consider the metric induced on the parameters by the scalar product on the field theory Hilbert space. Taking into account this metric, which is different from the naive Euclidean metric, one can then perform the optimization on RCMPS states viewed as a Riemannian manifold. Crucially, the induced metric is efficiently computable and efficiently invertible, and thus Riemannian techniques are not slower per iteration than their naive Euclidean counterparts. Ultimately, one can even use Riemannian quasi-Newton methods (implemented e.g. in OptimKit.jl) to further speed up convergence. The details of the numerical procedure can be found in [56].

## 3.5 Introducing defects

We now explain how to introduce extended operators (or defects) within the previous formalism. Let us first consider the expectation value of the extended operator $D_L = e^{-\mu \int_{-L}^{0} \hat{\phi}}$. We recall that while the defect has a finite size $L$, we took the thermodynamic limit for the bulk. Given the RCMPS ground state $|Q, R\rangle$ of the $\phi^4$ model, the expectation value of $D_L$ is given by:

$$ \langle 0, g | D_L | 0, g \rangle \approx \langle Q, R | e^{-\mu \int_{-L}^{0} \hat{\phi}} | Q, R \rangle. \qquad (3.16) $$

Here, $|Q, R\rangle \approx |0, g\rangle$ is the RCMPS approximation of a given bond dimension $D$ to the true ground state of the $\phi^4$ model at coupling $g$. Using the Baker-Campbell-Hausdorff formula
to normal order $\mathcal{D}_L$ (see appendix A.1) we have

$$ \langle Q, R | \mathcal{D}_L | Q, R \rangle = \mathcal{Z}_{-\mu G, -\mu G} \exp \left[ \frac{\mu^2}{2} \int_{\mathbb{R}} dx G(x)^2 \right], \quad (3.17) $$

with the modified “source”

$$ G(x) = \int_{-L}^{0} dy J(x-y), \quad (3.18) $$

and $J(x)$ being the source defined in eq. (3.8). As before, the path ordered exponential in $\mathcal{Z}_{-\mu G, -\mu G}$ is the solution to an ODE (3.15) now involving the modified source $G(x)$. The additional exponential factor in $\langle \mathcal{D}_L \rangle_{Q,R}$ can be incorporated into the ODE as a term proportional to the identity to get (in superoperator form)

$$ \frac{d}{dx} \rho(x) = \mathcal{L} \cdot \rho(x) - \mu G(x) [R\rho(x) + \rho(x)R^\dagger] + \frac{\mu^2}{2} G^2(x) \rho(x), \quad (3.19) $$

which gives

$$ \langle Q, R | \mathcal{D}_L | Q, R \rangle = \lim_{x \to \infty} \mathrm{tr}[\rho(x)], \quad \text{for the initial condition} \quad \lim_{x \to -\infty} \rho(x) = \rho_0. \quad (3.20) $$

Hence, computing the expectation value of the defect operator, which is extended, is not more difficult (asymptotically in $D$) than computing expectation values of local operators, as we advertised.

Correlation functions of local operators in the defect theory can be computed in a similar manner. As in section 3.3, we define the expectation value of a vertex operator in the full defect theory as

$$ \langle V_b(x) \rangle_{\text{defect}} = \frac{\langle 0, g | e^{b\phi(x)} : e^{-\mu \int_{-L}^{0} \hat{\phi}} | 0, g \rangle}{\langle 0, g | e^{-\mu \int_{-L}^{0} \hat{\phi}} | 0, g \rangle}, \quad (3.21) $$

where $|0, g\rangle$ is the exact ground state of the Hamiltonian at coupling $g$. This expression, upon iterated differentiation with respect to $b$, allows one to compute expectation value of $: \phi^k(x) :$ in the defect theory. We just computed the denominator, and the numerator can be evaluated in a similar way (see A.2.1 for more details). Ultimately, we obtain

$$ \frac{\langle Q, R | V_b(x) \mathcal{D}_L | Q, R \rangle}{\langle Q, R | \mathcal{D}_L | Q, R \rangle} = \frac{\mathcal{Z}_{s_x, s_x}}{\mathcal{Z}_{-\mu G, -\mu G}} \times \exp \left( -b\mu \int_{\mathbb{R}} dy \int_{-L}^{0} dz J(x-y)J(y-z) \right), \quad (3.22) $$

where

$$ s_x(y) = b J(x-y) - \mu \int_{-L}^{0} dz J(z-y). \quad (3.23) $$

Again, this means that the expectation value of the vertex operator in the defect theory can be computed efficiently as the ratio of the trace of the solution of matrix ODEs. Numerically, the only new difficulty is that r.h.s. of (3.22) is the ratio of two terms that grow exponentially in $L$, the defect size. This can however be alleviated with simple normalization strategies outlined in appendix A.2.3, which allow one to take the large defect limit $L \to \infty$ numerically.
# 4 Results

We can now numerically compute the expectation values of the defect operator $D_L$ and the field expectation value $\phi(x)$ with RCMPS for various values of the couplings $g$ and $\mu$. Since the physics only depends on $g/m^2$ and $\mu/m^2$, we fix $m = 1$ in all the RCMPS simulations, without loss of generality. First, we recall some perturbation theory basics that allow us to benchmark the RCMPS results in the weak coupling regime $g = 0.1$. Then, we move to the strong coupling regime, $g = 2$, where we can only compare RCMPS with itself. Subsequently, we increase the coupling to its critical value, $g_c \simeq 2.771525$ [60, 69], a regime that is particularly difficult with RCMPS at fixed bond dimension, but allows us to have some comparisons with conformal field theory results. Finally, we push the coupling to the symmetry-broken phase, $g = 4$, a regime where RCMPS allow us to consider two fundamentally different types of symmetry breaking: a natural one, following the defect perturbation, and a “frustrated” one, where the bulk is symmetry broken in the opposite way as the defect.

## 4.1 Perturbation theory calculations

Perturbation theory for the plain $\phi^4$ model (without the defect) is standard and explained in most textbooks [71, 72]. An advantage of the 1 + 1 dimensional case is that only the tadpole graphs are divergent, and they are removed by normal-ordering. Hence, all the Feynman diagrams appearing in the expansion are finite, and can be computed numerically e.g. with Monte Carlo [73].

To understand what happens when we add a defect, we first consider the free bulk case, i.e. $g = 0$. Since the defect is linear in the field, the whole functional integral is still Gaussian and one can compute it exactly. But, anticipating the interacting case, it still helps to see the result diagrammatically

For the defect expectation value in the free model, i.e. $g = 0$ and any $\mu$, we have that:

$$ \langle 0, 0 | D_L | 0, 0 \rangle = 1 + \sum_{k>0} \frac{(-\mu)^k}{k!} \langle 0, 0 | \left( \prod_{j=1}^{k} \int_{-L}^{0} dx_j \hat{\phi}(x_j) \right) | 0, 0 \rangle, \quad (4.1) $$

where $|0, 0\rangle$ is the ground state at $g = 0$, i.e. the free Fock vacuum of $H_m$. By systematically applying Wick's theorem, we can rewrite the above as a sum of products of integrated free massive propagators. Explicitly, to all orders in $\mu$ we have that:

$$ \begin{aligned} \langle 0, 0 | D_L | 0, 0 \rangle &= 1 + \frac{\mu^2}{2} \int_{-L}^{0} \int_{-L}^{0} dx_1 dx_2 G(x_1 - x_2; m) + \dots \\ &= \exp \left[ \frac{\mu^2}{2} \int_{-L}^{0} \int_{-L}^{0} dx_1 dx_2 G(x_1 - x_2; m) \right] \end{aligned} \quad (4.2) $$

with the free massive propagator given by

$$ G(x_i - x_j; m) := \langle 0, 0 | \hat{\phi}(x_i) \hat{\phi}(x_j) | 0, 0 \rangle = \frac{1}{2\pi} K_0(m|x_i - x_j|), \quad (4.3) $$
where $K_0$ is the modified Bessel function of the second kind. In terms of the following Feynman rules

$$ \begin{gathered} \overline{x_1 \quad x_2} \equiv G(x_1, x_2; m), \\ \multicolumn{1}{c}{\text{---}} \quad \text{---} \equiv -\mu \int_{-L}^{0} dx, \end{gathered} \qquad (4.4) $$

we can rewrite (4.2) as

$$ \begin{aligned} \langle 0, 0 | \mathcal{D}_L | 0, 0 \rangle &= \text{---} + \frac{1}{2} \text{---} + \frac{1}{8} \text{---} + \frac{1}{48} \text{---} + \dots \\ &= \exp \left[ \frac{1}{2} \text{---} \right]. \end{aligned} \qquad (4.5) $$

For the interacting model with $g \neq 0$ we proceed similarly, expanding both in $g$ and $\mu$. This time, the Feynman rules include bulk vertices

$$ \multicolumn{1}{c}{\text{---}} \quad \text{---} \equiv -4! g \int d^2 x, \qquad (4.6) $$

as well as propagators connecting the bulk and points on the defect. These propagators are given in eq. (4.3), where now each point can be in the bulk or on the defect.

At any fixed order in $g$, there are infinitely-many “disconnected” diagrams (i.e. where bulk and defect vertices are not connected to each other by any propagator), that “dress” the defect identity by defect interactions, just as in (4.5). At any fixed order in $g$, each connected diagram is “dressed” by the same factor, and there are finitely-many connected diagrams. In other words, and as far as connected correlators are considered, at any fixed order in $g$ the $\mu$-expansion truncates. Crucially, as we will later observe, this does not imply that the expansion is accurate at a fixed small $g$ and arbitrarily large $\mu$!

For the vacuum expectation value of $\mathcal{D}_L$ we find (neglecting $O(g^5)$ corrections):

$$ \begin{aligned} \frac{\langle 0, g | \mathcal{D}_L | 0, g \rangle}{\langle 0, 0 | \mathcal{D}_L | 0, 0 \rangle} &= 1 + \frac{1}{24} \text{---} + \frac{1}{12} \text{---} + \frac{1}{16} \text{---} + \frac{1}{72} \text{---} \\ &\quad + \frac{1}{24} \text{---} + \frac{1}{8} \text{---} + \frac{1}{32} \text{---} + \frac{1}{8} \text{---} \\ &\quad + \frac{1}{36} \text{---} + \frac{1}{24} \text{---} + \frac{1}{48} \text{---} + \frac{1}{144} \text{---} \\ &\quad + [34 \text{ diagrams } \mathcal{O}(g^4)]. \end{aligned} \qquad (4.7) $$

For the one-point function we find:

$$ \begin{aligned} \langle \phi(x) \rangle_{\text{defect}} \equiv \text{---} + \frac{1}{6} \text{---} + \frac{1}{6} \text{---} + \frac{1}{4} \text{---} + \frac{1}{12} \text{---} \\ + [13 \text{ diagrams } \mathcal{O}(g^3)] + [65 \text{ diagrams } \mathcal{O}(g^4)]. \end{aligned} \qquad (4.8) $$

We compute the diagrams in the expansions (4.7) and (4.8) up to order $g^4$ with an adaptive Monte Carlo method [73], at an accuracy sufficient to make the integration error (presumably) far lower than the neglected order $g^5$ terms.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Perturbative prediction for $\Gamma_0 \equiv \langle 0,g | \mathcal{D}_L | 0,g \rangle / \langle 0,0 | \mathcal{D}_L | 0,0 \rangle - 1$ up to $\mathcal{O}(g^4)$ vs. RCMPS (colored markers). Monte Carlo error bars are smaller than the size of the data points.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Up to $\mathcal{O}(g^4)$ perturbative prediction for $\equiv \langle\phi(x)\rangle_{\text{defect}} - \langle\phi(x)\rangle_{\text{defect},g=0}$ vs. RCMPS (colored markers). Monte Carlo error bars are smaller than the thickness of the lines.</figcaption></figure>

## 4.2 Weak coupling $g = 0.1$, comparison with perturbation theory

We first compare the perturbative calculation for the defect expectation value (4.7) and RCMPS in figure 3, where we have chosen $g = 0.1$ and considered two values $\mu = 1$ and $\mu = 5$ for the defect coupling. For $\mu = 1$, we observe that a modest bond dimension is sufficient to get well-converged RCMPS results, whereas we have to go to a fairly high order, $\mathcal{O}(g^4)$, in the perturbation theory to get similar values. For larger $\mu$, the defect is no longer a small perturbation as $L$ grows, and the perturbation theory results are already fairly inaccurate at $\mathcal{O}(g^4)$ (at least without Borel-Padé resummation) while RCMPS results remain just as good.

Next we consider the one-point function $\langle\phi(x)\rangle_{\text{defect}}$, from which we subtract the free theory contribution, in figure 4, for varying $x$, and fixed moderate defect size $L = 1$. We make similar observations. First, at small $\mu$, RCMPS quickly converges, while only the highest order in perturbation theory quantitatively agrees. Second, at large $\mu$, perturbation theory does not seem to converge, but RCMPS results converge even faster.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** (Rescaled logarithm of the) Defect vacuum expectation value from RCMPS computation as a function of the defect length $L$ in the strong-coupling, symmetric regime for $\mu = 1$ (left) and $\mu = 4$ (right).</figcaption></figure>

From the fact that the $g$-expansion is *exact* in $\mu$, which is a consequence of working with a defect that is linear in $\phi$, one might have expected perturbation theory to be accurate for fixed $g = 0.1$ and large $\mu$. However, as is clear from these numerical results, this is not the case. In fact, from the Feynman diagram expansion, it is clear that for very small $g$, the expansion is inaccurate when $\mu$ is large. For example, for the case of defect expectation values in (4.7), we see that diagrams of order $n$ scale like $O(g^n \mu^{2n+2})$, and thus become large when $\mu$ increases as $g$ is kept small but fixed.

## 4.3 Strong coupling $g = 2.0$ self-comparison

We now explore the strong-coupling (but still symmetric) regime $g = 2$. Using RCMPS, we can compute the vacuum expectation value of $D_L$ for finite (and not necessarily small) values of $g$ and $\mu$, and as a function of $L$. Since the defect expectation value grows exponentially as a function of $L$, we plot $L^{-1} \log\langle 0, g | D_L | 0, g \rangle$ in figure 5 for $\mu = 1$ and $\mu = 4$, which should go to a constant for large $L$.

As expected, the results converge quickly again, even for large $L$ and large $\mu$, in a regime where plain perturbation theory would be completely hopeless.

Next, we consider $|\langle\phi(x)\rangle_{\text{defect}}|$ for finite values of $g$, $\mu$ and $L$, as a function of $x$. For definiteness, in figure 6 we show the result corresponding to the choice $g = 2$, $\mu = 1, 4$, and $L = 6$. The results agree with our physical intuition. First, figure 6 shows that $|\langle\phi(x)\rangle_{\text{defect}}|$ decays exponentially as we move away from the center of the defect $x = -L/2$ into the bulk, as expected. Second, increasing the defect coupling $\mu$ (keeping all other parameters fixed), $|\langle\phi(x)\rangle_{\text{defect}}|$ increases, as it should since the bulk field tends to align with the defect. Third, the field expectation value remains continuous at $x = 0$, i.e. as we cross the defect. For $x < 0$ it measures a *defect magnetization* and stabilizes to a constant value in the middle of the defect.

At large $\mu$, the RCMPS results, although accurate at $D \ge 24$, are arguably slower to converge, especially inside the defect. We believe this is fairly intuitive, and a property of all variational methods. Our RCMPS states are approximate ground states of the bulk
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Bulk one-point functions from RCMPS computation as a function of the distance from the defect (shaded region) in the strong-coupling, symmetric regime, for $\mu = 1$ (left) and $\mu = 4$ (right).</figcaption></figure>

model, obtained by minimizing an energy density. From the ground state, one can define a probability distribution for the spatial average of the field in any finite size interval. The RCMPS approximates this probability very well overall, but is inevitably less accurate in the tails, for very small probabilities that almost do not contribute to the energy density. As we increase $\mu$ we shift the center of the distribution, and thus get more sensitive to its tail. As a result, we need a larger bond dimension to preserve the accuracy. This shift of the distribution is stronger closer to the defect, and even more so inside it, and thus we need larger bond dimensions inside the defect. Note however that even deep in the defect and at large $\mu$ we can easily reach bond dimensions such that the results are visually converged in figure 6. We do not know of any other method that could reach a similar accuracy for such large perturbations both in $g$ and $\mu$ in the continuum.

## 4.4 Critical coupling $g \simeq 2.771525$, comparison with scaling expectations

We now tune the bulk model to its critical coupling, $g_c \simeq 2.771525$ [60, 69]. In this strong coupling regime, one might be able to test RCMPS not just against itself, but also against universal scaling predictions from Conformal Field Theory. Indeed, at criticality, the $\phi^4$ model is in the Ising universality class, which is well understood.

However, this universal regime is also where the accuracy of RCMPS degrades because of their finite entanglement entropy. In principle, one could improve the accuracy of the RCMPS results by doing a finite entanglement scaling of the results [59], something we keep for future work. Here, we stick with a plain variational study, which also helps to demonstrate the limitations of RCMPS.

We first consider the expectation value of the defect. From conformal field theory, the expectation value of a defect line of length $L$ in a critical theory should scale like the two-point function of a *defect creation operator*, inserted at the line defect endpoints (see e.g. the appendix of ref. [74]). In our case, if we identify $\phi$ with the spin $\sigma$ in the 2d Ising model,
<figure><img src="image_9.png" /><figcaption>Figure 7. Defect vacuum expectation value from RCMPS computation as a function of the defect length L at critical bulk coupling, for $\mu = 1$ and $\mu = 4$. The dashed line is the exact result given in eq. (4.10), shifted along the vertical axis by an (arbitrary) constant for visualization purposes.</figcaption></figure>

the corresponding scaling exponent $\Delta_{+0}$ was computed in [27] and we expect

$$ \langle 0, g_c | \mathcal{D}_L | 0, g_c \rangle \sim e^{-EL} L^{-2\Delta_{+0}}, \qquad (4.9) $$

when $L \to +\infty$ with $\Delta_{+0} = 1/32$ and E a non-universal term. Thus should have

$$ \left(1 - L \frac{d}{dL}\right) \log \langle 0, g_c | \mathcal{D}_L | 0, g_c \rangle \sim -\frac{1}{16} \log L. \qquad (4.10) $$

In figure 7, we compare this candidate scaling against RCMPS calculations, for $\mu = 1$ and $\mu = 4$. For sufficiently small L the RCMPS results are accurate and well converged, but convergence becomes dramatically slower at large distance as expected. Even pushing the bond dimension to $D = 64$, we do not see the expected large distance scaling very clearly.

It is possible that our bond dimension is just too low, or that we would need a refined entanglement scaling analysis. It is also possible that our scaling guess (4.9) is incorrect because of additional non-universal UV effects, or substantially corrected by less relevant operators in the IR which are not decayed enough for the lengths L we can reliably probe.

In figure 8 we show results for $|\langle\phi(x)\rangle_{\text{defect}}|$ with $\mu = 1, 4$ and $L = 30$. For $1 \ll x \ll L$ we should again be sensitive to the critical (scaling) behavior of the model. The scaling behavior of the one-point function was computed in [13] using Boundary Conformal Field Theory techniques (see [4, 9] for a review):

$$ \langle \phi(x) \rangle_{\text{defect}} = \frac{2^{1/8} \sqrt{C_{\phi}}}{|x|^{1/8}}, \qquad (4.11) $$

where $C_{\phi}$ is a non-universal normalization constant which is obtained from the two-point function of the field (without defect) at large distance

$$ \langle 0, g_c | \hat{\phi}(x_1) \hat{\phi}(x_2) | 0, g_c \rangle \sim \frac{C_{\phi}}{|x_1 - x_2|^{1/4}}. \qquad (4.12) $$
<figure><img src="image_10.png" /><figcaption>**Figure 8.** One-point functions obtained from RCMPS computation as a function of the distance from the defect at critical coupling, for $\mu = 1$ and $\mu = 4$. The dashed line is the expected exact scaling given in eq. (4.11).</figcaption></figure>

We estimate it crudely in appendix B by fitting the asymptotic behavior of the two-point function without the defect and get $C_{\phi} \simeq 0.34$. In figure 8 we observe a slightly better qualitative agreement with the scaling prediction (4.12) than in the pure defect case considered before. However, without a careful finite entanglement scaling, our bond dimension remains too low to reliably fit an exponent.

## 4.5 Symmetry broken phase, coupling $g = 4$

We now go deep in the symmetry broken phase. There are now two options for symmetry breaking, that correspond to very different physical situations:

* The bulk can be symmetry broken in the direction favored by the defect, $\text{sign}(\langle 0, g | \hat{\phi}(x) | 0, g \rangle) = -\text{sign}(\mu)$. This is what would happen spontaneously in a statistical mechanical system with a defect, initialized in the symmetric phase but then slowly cooled down. The bulk would then pick the average field value favored by the defect.

* The bulk can be symmetry broken “against” the defect, i.e. $\text{sign}(\langle 0, g | \hat{\phi}(x) | 0, g \rangle) = +\text{sign}(\mu)$. Physically, this would occur if symmetry breaking happens first, yielding positive $\langle 0, g | \hat{\phi}(x) | 0, g \rangle$, and only then the defect coupling $\mu$ is progressively sent to positive values. The bulk field would then remain stuck “against” the defect.

In figure 9, we show the RCMPS predictions for $L^{-1} \log \langle D_L \rangle$ as a function of $L$ in these two cases. When the symmetry is broken along the defect, the results are fairly expected: the defect expectation value increases more or less exponentially with system size, i.e. $L^{-1} \log \langle D_L \rangle$ is almost constant and positive. RCMPS results are well converged already for very small bond dimensions.

However, when the defect is against the bulk field, the defect expectation value first decreases as a function of $L$ before increasing again, with an asymptotic exponential growth at large $L$. Convergence as a function of the bond dimension is slower, as expected from a
<figure><img src="image_11.png" /><figcaption>Figure 9. Defect vacuum expectation value from RCMPS computation as a function of the defect length L in the symmetry-broken regime. In (a) symmetry breaking along the defect, in (b) symmetry breaking against the defect. In both plots $\langle Q, R|\hat{\phi}(x)|Q, R\rangle > 0$.</figcaption></figure>

variational method: we probe field values that are very atypical compared to those in the (approximate) ground state. Nonetheless, for $D \le 24$ the results appear well converged.

This subtle behavior as a function of L in the second case is physically interesting, and can be understood with a bit of statistical field theory intuition. When the defect is small, every point inside the defect has a field environment dominated by the bulk. This is very unfavorable, because the defect penalizes precisely the field sign favored in the bulk. However, when the defect is sufficiently large, most points inside the defect have a field environment dominated by the defect itself, with the right sign, and thus contribute positively to the defect expectation value. Intuitively, this change of behavior can happen only if $\mu$ is large enough, otherwise the pull of the bulk always dominates.

To support this intuition, we consider the asymptotic growth rate of the defect expectation value:

$$ \lim_{L \to \infty} \frac{\log \langle 0, g | \mathcal{D}_L | 0, g \rangle}{L}, \qquad (4.13) $$

in the symmetry broken phase, as a function of $\mu$. Numerically, we approximate this $L \to \infty$ limit by $L = 32$, which is much larger than the inverse mass scale. The corresponding results are shown in figure 10.

We observe that the asymptotic growth rate of the defect expectation value is an approximate parabola that is not centered on zero but on positive $\mu$, as a clear sign of symmetry breaking. It is positive for all negative $\mu$ (this corresponds to symmetry breaking “along” the defect) as expected, negative for small positive $\mu$ (symmetry breaking “against” the defect, and a defect environment dominated by the bulk influence), and finally positive when $\mu$ is large enough (symmetry breaking “against” the defect, and a defect environment dominated by the defect influence).

We finally consider the 1-point function in the presence of a defect of fixed size $L = 6$ in figure 11, where we used again RMCPS states with $\langle Q, R|\hat{\phi}(x)|Q, R\rangle > 0$. These results
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Asymptotic growth rate of the defect expectation value as a function of $\mu$. The results are shown for RCMPS with positive field expectation value $\langle Q, R|\hat{\phi}(x)|Q, R\rangle > 0$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** One-point functions from RCMPS computation as a function of the distance from the defect (shaded region) in the symmetry-broken regime. In (a) symmetry breaking along the defect, in (b) symmetry breaking against the defect. In both plots $\langle Q, R|\hat{\phi}(x)|Q, R\rangle > 0$.</figcaption></figure>

show the expected sign change in the field when moving from the defect to the bulk when the bulk is against the defect.

## 5 Conclusion

In this article, we accurately computed the expectation values in a 1 + 1 dimensional scalar field theory with a magnetic line defect, using a particular variational method (RCMPS). The method is impressively efficient in many cases, albeit with a few limitations.

First, and most importantly, we could carry non-perturbative calculations without any discretizaion or needing a UV cutoff. We could also easily take the thermodynamic limit in the bulk first, for a defect of fixed size, without extra numerical cost.

Second, we could consider a wide range of bulk couplings $g$, and go deep in the strong coupling and even symmetry-broken regime without any changes to our ansatz. We could
also consider large defect coupling $\mu$, which as we recall, is impossible in perturbation theory even at fixed $g \ll 1$. Taking large defects, numerically $L$ much larger than other length-scales, brought essentially no drop in precision, and only a marginal increase in computational cost.

Putting aside the UV problems, some of the expectation values we have computed would have been difficult to estimate even on a finite lattice. Indeed, the symmetry-broken regime where the defect pulls “against” the bulk average, would be non-trivial to probe with the Monte-Carlo method. First, because one would have to take limits carefully to ensure the bulk remains broken against the defect during sampling. Second, because the defect expectation values are dominated by extremely rare events in that regime.

Our method shines comparatively less to estimate universal properties, where mastering the fine-grained details of the continuous model is not needed. This is particularly the case at the critical point, where our entanglement cutoff introduces an effective length-scales that breaks the scale invariance of the correlators. To estimate universal exponents, one is certainly better off with a (not necessarily fine) lattice discretization, solved with standard MPS. Universality will guarantee that the infrared physics remains identical, but the much lower computational cost of MPS would allow to reach far larger bond dimensions, and thus obtain much better precision in the IR.

There are natural continuations of the present work. First, one could extend finite entanglement scaling techniques developed on the lattice to the continuum. This would certainly provide much better estimates of universal and non-universal quantities at criticality (even if, for universal quantities, lattice discretization will likely still dominate). Second, one could consider magnetic defects with imaginary coupling. While this presents no difficulty with RCMPS (one can take $\mu$ imaginary in all the formulas we presented), the lattice Monte-Carlo method is effectively unusable because of a growing sign problem as a function of $L$. Considering other defects (non-magnetic) like $\phi^2$ would be interesting as well, but we have not found a numerically efficient method with RCMPS yet. Finally, we could apply the present method to other bosonic field theories (like the $\phi^6$, Sine-Gordon, and Sinh-Gordon models) or even to models with multiple scalar fields (like the $O(2)$ model).

## Acknowledgments

We thank G. Cuomo for discussions. We are supported by the European Union (ERC, QFT.zip project, Grant Agreement no. 101040260). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.

## A Derivation of RCMPS formulae for defects

In this appendix, we present the full derivation of the RCMPS formulae given in section 3.5. Below $|Q, R\rangle$ is a RCMPS approximation to the true (interacting) ground state $|0, g\rangle$ of the $\phi^4$ model at coupling $g$.
## A.1 Expectation value of the defect operator

We start with the expectation value of the defect operator $\mathcal{D}_L$ defined in (2.7)

$$ \langle 0, g | \mathcal{D}_L | 0, g \rangle = \frac{\int [D\phi] \mathcal{D}_L e^{-S_B}}{\int [D\phi] e^{-S_B}}. \qquad (A.1) $$

Plugging in the definition of $\hat{\phi}(x)$

$$ \hat{\phi}(x) = \int_{\mathbb{R}} dy J(x-y)[\hat{a}(y) + \hat{a}^{\dagger}(y)], \quad J(x) = \frac{1}{2\pi} \int_{\mathbb{R}} \frac{dk}{\sqrt{2\omega_k}} e^{-ikx}. \qquad (A.2) $$

we can use the Baker-Campbell-Hausdorff (BCH) formula to put $\mathcal{D}_L$ in normal-ordered form

$$ \begin{aligned} \mathcal{D}_L &=: \mathcal{D}_L : \times \exp \left[ \frac{\mu^2}{2} \int_{[-L,0]^2} dx dx' \int_{\mathbb{R}} dy J(x-y) J(x' - y) \right] \\ &=: \mathcal{D}_L : \times \langle 0, 0 | \mathcal{D}_L | 0, 0 \rangle, \end{aligned} \qquad (A.3) $$

where $|0, 0\rangle$ is the ground state at coupling 0, hence the Fock vacuum.

So far, everything was exact. Replacing $|0, g\rangle$ by its RCMPS approximation we get

$$ \langle 0, g | : \mathcal{D}_L : | 0, g \rangle \approx \langle Q, R | : \mathcal{D}_L : | Q, R \rangle = \mathcal{Z}_{-\mu G, -\mu G}, \qquad (A.4) $$

and thus (A.3) gives:

$$ \langle 0, g | \mathcal{D}_L | 0, g \rangle \approx \mathcal{Z}_{-\mu G, -\mu G} \times \exp \left[ \frac{\mu^2}{2} \int_{\mathbb{R}} dy G(y)^2 \right], \qquad (A.5) $$

where $\mathcal{Z}_{j',j}$ is the RCMPS generating functional defined in (3.7), and $G$ is the modified “source”, i.e.

$$ \mathcal{Z}_{-\mu G, -\mu G} = \mathrm{tr} \left[ \mathcal{P} \exp \int_{\mathbb{R}} dx \, \mathbb{T} - \mu G(x) (R \otimes 1 + 1 \otimes R^*) \right] \quad \text{with} \quad G(x) := \int_{-L}^{0} dy J(x-y). \qquad (A.6) $$

This path ordered exponential can be written as the solution to an ordinary differential equation (ODE). Putting the ODE in superoperator form as in the main text, we finally get

$$ \langle 0, g | \mathcal{D}_L | 0, g \rangle \approx \lim_{x \to \infty} \mathrm{tr} [\rho(x)], \quad \text{for the initial condition} \quad \lim_{x \to -\infty} \rho(x) = \rho_0. \qquad (A.7) $$

and

$$ \frac{d}{dx} \rho(x) = \mathcal{L} \cdot \rho(x) - \mu G(x) [R\rho(x) + \rho(x)R^{\dagger}] + \frac{\mu^2}{2} G^2(x) \rho(x), \qquad (A.8) $$

where we recall that

$$ \mathcal{L} \cdot \rho = Q \rho + \rho Q^{\dagger} + R \rho R^{\dagger}. \qquad (A.9) $$
## A.2 One-point functions in the defect model

### A.2.1 Vertex operators in the defect model

Next, we consider the vertex operator in the full defect theory (3.21), which we differentiate in A.2.2 to obtain expectation values of one point functions in the $\phi^4$ model in presence of the defect. It is defined as

$$ \langle V_b(x) \rangle_{\text{defect}} = \frac{\langle 0, g | e^{b\hat{\phi}(x)} : e^{-\mu \int_{-L}^{0} \hat{\phi}} | 0, g \rangle}{\langle 0, g | e^{-\mu \int_{-L}^{0} \hat{\phi}} | 0, g \rangle}. \qquad (A.10) $$

We first focus on the numerator, as the denominator has just been computed.

The first step, again exact, is to normal-order the operator being evaluated using the Baker-Campbell-Hausdorff formula:

$$ : e^{b \hat{\phi}(x)} : e^{-\mu \int_{-L}^{0} dy \, \hat{\phi}(y)} = : e^{b \hat{\phi}(x) - \mu \int_{-L}^{0} dy \, \hat{\phi}(y)} : \langle 0, 0 | \mathcal{D}_L | 0, 0 \rangle \\ \times \exp \left[ -b \mu \int_{\mathbb{R}} dy \, J(x-y) G(y) \right], \qquad (A.11) $$

where we have explicitly separated the $\langle 0, 0 | \mathcal{D}_L | 0, 0 \rangle$ for convenience as it also appears in the denominator and thus ultimately cancels out.

Replacing the true ground state by its variational approximation yields

$$ \langle 0, g | : e^{b\hat{\phi}(x)} : e^{-\mu \int_{-L}^{0} \hat{\phi}} | 0, g \rangle \approx \mathcal{Z}_{s_x, s_x} \langle 0, 0 | \mathcal{D}_L | 0, 0 \rangle \\ \times \exp \left[ -b \mu \int_{\mathbb{R}} dy \, J(x-y) G(y) \right], \qquad (A.12) $$

where

$$ s_x(y) := b J(x-y) - \mu \int_{-L}^{0} dz \, J(z-y) = b J(x-y) - \mu G(y). \qquad (A.13) $$

Again, we can write $\mathcal{Z}_{s_x, s_x}$ as the solution of an ODE to get

$$ \langle Q, R | : e^{b \hat{\phi}(x)} : \mathcal{D}_L | Q, R \rangle = \langle 0, 0 | \mathcal{D}_L | 0, 0 \rangle \lim_{y \to +\infty} \text{tr}[\rho(y)], \qquad (A.14) $$

with

$$ \frac{d}{dy} \rho(y) = \mathcal{L} \cdot \rho(x) + s_x(y) [R\rho(x) + \rho(x)R^\dagger] - b\mu G(y)J(x-y)\rho(y), \qquad (A.15) $$

and with the initial condition $\lim_{y \to -\infty} \rho(y) = \rho_0$ with $\text{tr}[\rho_0] = 1$. Finally, we can put numerator and denominator together, which cancels the $\langle 0, 0 | \mathcal{D}_L | 0, 0 \rangle$ term and we get

$$ \langle V_b(x) \rangle_{\text{defect}} \approx \lim_{y \to \infty} \frac{\text{tr}[\rho(y)]}{\text{tr}[\tilde{\rho}(y)]}, \qquad (A.16) $$

where $\tilde{\rho}(y)$ is the solution of (A.15) with $b = 0$ (or equivalently of (A.8) without $G^2$ term).
### A.2.2 Field monomials in the defect model

With explicit expressions for defect vertex operators, we may now compute expectation values of field monomials in the full defect theory. The latter are defined as:

$$ \langle :\phi^n(x): \rangle_{\text{defect}} = \frac{\langle 0, g | :\hat{\phi}^n(x) : e^{-\mu} \int_{-L}^{0} dx \, \hat{\phi}(x) | 0, g \rangle}{\langle 0, g | e^{-\mu} \int_{-L}^{0} dx \, \hat{\phi}(x) | 0, g \rangle} \quad (A.17) $$

Again, we already have the denominator from (A.1). For the numerator, we can differentiate the defect vertex operator w.r.t. $b$,

$$ \langle 0, g | : \hat{\phi}^n(x) : e^{-\mu} \int_{-L}^{0} dx \, \hat{\phi}(x) | 0, g \rangle = \frac{\partial^n}{\partial b^n} \langle 0, g | : e^{b} \hat{\phi}(x) : \mathcal{D}_{L} | 0, g \rangle \Big|_{b=0} \quad (A.18) $$

The strategy is now to forward differentiate the ordinary differential equation giving us $\rho(y)$ appearing in (A.16)

$$ \langle 0, g | : \hat{\phi}^n(x) : e^{-\mu} \int_{-L}^{0} dx \, \hat{\phi}(x) | 0, g \rangle \approx \langle Q, R | : \hat{\phi}^n(x) : e^{-\mu} \int_{-L}^{0} dx \, \hat{\phi}(x) | Q, R \rangle \quad (A.19) $$

$$ = \langle 0, 0 | \mathcal{D}_L | 0, 0 \rangle \times \lim_{y \to +\infty} \text{tr} \left[ \rho^{(n)}(y) \right], \quad (A.20) $$

where $\rho^{(n)}(y) := \partial_b^n \rho(y)|_{b=0}$ and $\rho$ obeys the ODE (A.15). The matrices $\rho^{(k)}(y)$ for $0 \le k \le n$ obey the triangular system of matrix ODEs

$$ \frac{d}{dy} \rho^{(k)}(y) = \mathcal{L} \cdot \rho^{(k)}(y) + kJ(x-y) \left[ R\rho^{(k-1)}(y) + \rho^{(k-1)}(y)R^{\dagger} \right] \quad (A.21) $$

$$ - \mu G(y) \left[ R \rho^{(k)}(y) + \rho^{(k)}(y) R^{\dagger} \right] - k \mu J(x-y) G(y) \rho^{(k-1)}(y), \quad (A.22) $$

with the initial conditions, $\lim_{y \to -\infty} \rho^{(0)}(y) = \rho_0$ and $\lim_{y \to -\infty} \rho^{(n)}(y) = 0$ for $k > 0$. Here, again the factor $\langle 0, 0 | \mathcal{D}_L | 0, 0 \rangle$ arises both in the numerator (A.20) and the denominator (A.3) and hence drops out in the computation of monomials (A.17). Putting all together, we find

$$ \langle :\phi^n(x): \rangle_{\text{defect}} = \lim_{y \to +\infty} \frac{\text{tr}[\rho^{(n)}(y)]}{\text{tr}[\rho^{(0)}(y)]} \quad (A.23) $$

### A.2.3 Normalization tricks

The vertex operator and field monomial expectation values in the presence of the defect given in (A.16) and (A.23) are ratios of traces, each of which is typically growing exponentially with the size $L$ of the defect. Thus, computing each term independently and taking the ratio quickly gives large numerical errors in double precision, especially when numerically taking the semi-infinite defect line limit. This should be avoided.

A simple way to make things better behaved is to periodically normalize both the matrix in the numerator and the matrix in the denominator by the trace of the denominator as the ODE is being solved. For example, in the case of the computation of the field expectation

$$ \langle :\phi(x): \rangle_{\text{defect}} = \lim_{y \to +\infty} \frac{\text{tr}[\rho^{(1)}(y)]}{\text{tr}[\rho^{(0)}(y)]}, \quad (A.24) $$
we can pick a number of checkpoints $y_k$ (either a priori, or each time $\text{tr}[\rho^{(0)}(y)]$ gets too large or too small) and renormalize both $\rho^{(0)}(y_k)$ and $\rho^{(1)}(y_k)$:

$$ \rho^{(0)}(y_k^+) = \frac{\rho^{(0)}(y_k)}{\text{tr}[\rho^{(0)}(y_k)]}, \qquad (A.25) $$

$$ \rho^{(1)}(y_k^+) = \frac{\rho^{(1)}(y_k)}{\text{tr}[\rho^{(0)}(y_k)]}. \qquad (A.26) $$

Inserting these normalizations do not change the value of the expectation value (A.24).

Alternatively, one could define the continuously normalized matrices

$$ \sigma^{(0)}(y) = \frac{\rho^{(0)}(y)}{\text{tr}[\rho^{(0)}(y)]} \qquad (A.27) $$

$$ \sigma^{(1)}(y) = \frac{\rho^{(1)}(y)}{\text{tr}[\rho^{(0)}(y)]}, \qquad (A.28) $$

which indeed obey a system of closed ODEs. However these ODEs are non-linear, and we found it numerically preferable to use the previous periodic normalization strategy. No matter which strategy is used, one can compute expectation values for large or even numerically semi-infinite defects without suffering from instabilities of precision loss.

### A.3 Strategy for more general expectation values

Following the same techniques, one can compute expectation values of arbitrary strings of vertex operators or normal-ordered monomials in the presence of the defect. The steps are the same as before, and we only outline the general strategy here. One first considers expectation values of products of vertex operators

$$ \langle V_{b_1}(x_1) V_{b_2}(x_2) \dots V_{b_n}(x_n) \rangle_{\text{defect}} = \frac{\langle 0, g | e^{b_1 \hat{\phi}(x_1)} : e^{b_2 \hat{\phi}(x_2)} : \dots : e^{b_n \hat{\phi}(x_n)} : e^{-\mu \int_{-L}^{0} \hat{\phi}} | 0, g \rangle}{\langle 0, g | e^{-\mu \int_{-L}^{0} \hat{\phi}} | 0, g \rangle}, \qquad (A.29) $$

for $x_1 < x_2 < \dots < x_n$. Then, focusing on the numerator, we use the Baker-Campbell-Hausdorff formula iteratively to put all the fields under the same exponential

$$ \begin{gathered} : e^{b_1 \hat{\phi}(x_1)} : : e^{b_2 \hat{\phi}(x_2)} : \dots : e^{b_n \hat{\phi}(x_n)} : e^{-\mu \int_{-L}^{0} \hat{\phi}} = \\ \mathcal{N}(b_1, b_2, \dots, b_n, x_1, x_2, \dots, x_n) \times : e^{b_1 \hat{\phi}(x_1) + b_2 \hat{\phi}(x_2) + \dots + b_n \hat{\phi}(x_n)} : e^{-\mu \int_{-L}^{0} \hat{\phi}}, \end{gathered} \qquad (A.30) $$

where $\mathcal{N}$ is an explicit scalar contribution coming from the normal ordering. The expectation value of the remaining normal-ordered exponential on a RCMPS is then just $\mathcal{Z}_{s,s}$ with

$$ s(y) := \sum_{j=1}^{n} b_j J(x_j - y) - \mu G(y), \qquad (A.31) $$

and can thus be obtained by solving a simple linear ODE like before.
<figure><img src="image_14.png" /><figcaption>Figure 12. In (a), up to $O(g^4)$ perturbative prediction for $\Gamma_2(x) \equiv \langle 0,0|\phi(0)\phi(x)|0,g \rangle - \langle 0,g|\phi(0)\phi(x)|0,0 \rangle$ vs. RCMPS (colored markers). Monte Carlo error bars are smaller than the size of the data points. In (b), no-defect two-point functions from RCMPS computation as a function of the distance from the defect at critical coupling. The dashed line is eq. (B.3) with $C_\phi \simeq 0.34$.</figcaption></figure>

## B.2 Critical coupling regime

We now tune the bulk coupling to its critical value, $g_c \simeq 2.771525$ [54–56, 69]. In figure 12 we show results for $\langle 0, g_c|\hat{\phi}(x)\hat{\phi}(0)|0, g_c \rangle$. In the scaling region the result is compatible with the power-law decay for critical Ising model [4], i.e.

$$ \langle \phi(x) \phi(0) \rangle_{\text{bulk}} = \frac{C_{\phi}}{|x|^{1/4}}, \qquad (B.3) $$

where $C_\phi \simeq 0.34$ is a non-universal coefficient which we evaluate numerically from the scaling region of the plot in figure 12.

**Data Availability Statement.** This article has no associated data or the data will not be deposited.

**Code Availability Statement.** This article has no associated code or the code will not be deposited.

**Open Access.** This article is distributed under the terms of the Creative Commons Attribution License (CC-BY4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.

## References

[1] K.G. Wilson, *The Renormalization Group: Critical Phenomena and the Kondo Problem*, Rev. Mod. Phys. 47 (1975) 773 [INSPIRE].

[2] J. Kondo, *Theory of Dilute Magnetic Alloys*, in *Solid State Physics*, F. Seitz, D. Turnbull and H. Ehrenreich eds., Elsevier (1970), p. 183–281 [DOI:10.1016/s0081-1947(08)60616-5].
[3] J. Cardy, *Scaling and Renormalization in Statistical Physics*, Cambridge University Press (1996) [DOI:10.1017/cbo9781316036440].

[4] P. Di Francesco, P. Mathieu and D. Senechal, *Conformal Field Theory*, Springer-Verlag, New York (1997) [DOI:10.1007/978-1-4612-2256-9] [INSPIRE].

[5] D. Gaiotto, A. Kapustin, N. Seiberg and B. Willett, *Generalized Global Symmetries*, JHEP 02 (2015) 172 [arXiv:1412.5148] [INSPIRE].

[6] C. Cordova, T.T. Dumitrescu, K. Intriligator and S.-H. Shao, *Snowmass White Paper: Generalized Symmetries in Quantum Field Theory and Beyond*, in the proceedings of the *Snowmass 2021*, Seattle, U.S.A., July 17-26 (2022) [arXiv:2205.09545] [INSPIRE].

[7] D.C. Lewellen, *Sewing constraints for conformal field theories on surfaces with boundaries*, Nucl. Phys. B 372 (1992) 654 [INSPIRE].

[8] I. Runkel, *Boundary structure constants for the A series Virasoro minimal models*, Nucl. Phys. B 549 (1999) 563 [hep-th/9811178] [INSPIRE].

[9] J.L. Cardy, *Boundary conformal field theory*, hep-th/0411189 [INSPIRE].

[10] P. Liendo, L. Rastelli and B.C. van Rees, *The Bootstrap Program for Boundary CFT$_d$*, JHEP 07 (2013) 113 [arXiv:1210.4258] [INSPIRE].

[11] M. Billò, V. Gonçalves, E. Lauria and M. Meineri, *Defects in conformal field theory*, JHEP 04 (2016) 091 [arXiv:1601.02883] [INSPIRE].

[12] E. Lauria, P. Liendo, B.C. Van Rees and X. Zhao, *Line and surface defects for the free scalar field*, JHEP 01 (2021) 060 [arXiv:2005.02413] [INSPIRE].

[13] G. Cuomo, Z. Komargodski and M. Mezei, *Localized magnetic field in the $O(N)$ model*, JHEP 02 (2022) 134 [arXiv:2112.10634] [INSPIRE].

[14] T. Nishioka, Y. Okuyama and S. Shimamori, *The $\epsilon$ expansion of the $O(N)$ model with line defect from conformal field theory*, JHEP 03 (2023) 203 [arXiv:2212.04076] [INSPIRE].

[15] A. Gimenez-Grau, *Probing magnetic line defects with two-point functions*, arXiv:2212.02520 [INSPIRE].

[16] S. Giombi, E. Helfenberger and H. Khanchandani, *Line defects in fermionic CFTs*, JHEP 08 (2023) 224 [arXiv:2211.11073] [INSPIRE].

[17] D. Rodriguez-Gomez and J.G. Russo, *Defects in scalar field theories, RG flows and dimensional disentangling*, JHEP 11 (2022) 167 [arXiv:2209.00663] [INSPIRE].

[18] L. Bianchi, D. Bonomi and E. de Sabbata, *Analytic bootstrap for the localized magnetic field*, JHEP 04 (2023) 069 [arXiv:2212.02524] [INSPIRE].

[19] O. Aharony et al., *Phases of Wilson Lines in Conformal Field Theories*, Phys. Rev. Lett. 130 (2023) 151601 [arXiv:2211.11775] [INSPIRE].

[20] G. Cuomo, Z. Komargodski, M. Mezei and A. Raviv-Moshe, *Spin impurities, Wilson lines and semiclassical*, JHEP 06 (2022) 112 [arXiv:2202.00040] [INSPIRE].

[21] O. Aharony et al., *Phases of Wilson lines: conformality and screening*, JHEP 12 (2023) 183 [arXiv:2310.00045] [INSPIRE].

[22] J. Barrat, P. Liendo and P. van Vliet, *Line defect correlators in fermionic CFTs*, arXiv:2304.13588 [INSPIRE].

[23] L. Bianchi, D. Bonomi, E. de Sabbata and A. Gimenez-Grau, *Analytic bootstrap for magnetic impurities*, JHEP 05 (2024) 080 [arXiv:2312.05221] [INSPIRE].
[24] I. Nagar, A. Sever and D.-L. Zhong, *Planar RG flows on line defects*, JHEP 06 (2024) 110
[arXiv:2404.07290] [INSPIRE].

[25] O. Diatlyk, H. Khanchandani, F.K. Popov and Y. Wang, *Defect fusion and Casimir energy in higher dimensions*, JHEP 09 (2024) 006 [arXiv:2404.05815] [INSPIRE].

[26] P. Kravchuk, A. Radcliffe and R. Sinha, *Effective theory for fusion of conformal defects*, arXiv:2406.04561 [INSPIRE].

[27] G. Cuomo, Y.-C. He and Z. Komargodski, *Impurities with a cusp: general theory and 3d Ising*, JHEP 11 (2024) 061 [arXiv:2406.10186] [INSPIRE].

[28] M. Vojta, C. Buragohain and S. Sachdev, *Quantum impurity dynamics in two-dimensional antiferromagnets and superconductors*, Phys. Rev. B 61 (2000) 15152.

[29] S. Sachdev and M. Vojta, *Quantum impurity in an antiferromagnet: Nonlinear sigma model theory*, Phys. Rev. B 68 (2003) 064419 [cond-mat/0303001] [INSPIRE].

[30] S. Florens, L. Fritz and M. Vojta, *Kondo Effect in Bosonic Spin Liquids*, Phys. Rev. Lett. 96 (2006) 036601.

[31] F.F. Assaad and I.F. Herbut, *Pinning the order: the nature of quantum criticality in the Hubbard model on honeycomb lattice*, Phys. Rev. X 3 (2013) 031010 [arXiv:1304.6340] [INSPIRE].

[32] F. Parisen Toldin, F.F. Assaad and S. Wessel, *Critical behavior in the presence of an order-parameter pinning field*, Phys. Rev. B 95 (2017) 014401 [arXiv:1607.04270] [INSPIRE].

[33] A. Allais, *Magnetic defect line in a critical Ising bath*, arXiv:1412.3449.

[34] A. Franchi, D. Rossini and E. Vicari, *Critical crossover phenomena driven by symmetry-breaking defects at quantum transitions*, Phys. Rev. E 105 (2022) 034139 [arXiv:2201.00390] [INSPIRE].

[35] D. Rodriguez-Gomez, *A scaling limit for line and surface defects*, JHEP 06 (2022) 071 [arXiv:2202.03471] [INSPIRE].

[36] F.K. Popov and Y. Wang, *Non-perturbative defects in tensor models from melonic trees*, JHEP 11 (2022) 057 [arXiv:2206.14206] [INSPIRE].

[37] A. Gimenez-Grau, E. Lauria, P. Liendo and P. van Vliet, *Bootstrapping line defects with $O(2)$ global symmetry*, JHEP 11 (2022) 018 [arXiv:2208.11715] [INSPIRE].

[38] W.H. Pannell and A. Stergiou, *Line defect RG flows in the $\varepsilon$ expansion*, JHEP 06 (2023) 186 [arXiv:2302.14069] [INSPIRE].

[39] L. Hu, Y.-C. He and W. Zhu, *Solving conformal defects in 3D conformal field theory using fuzzy sphere regularization*, Nature Commun. 15 (2024) 3659 [arXiv:2308.01903] [INSPIRE].

[40] Z. Zhou, D. Gaiotto, Y.-C. He and Y. Zou, *The $g$-function and defect changing operators from wavefunction overlap on a fuzzy sphere*, SciPost Phys. 17 (2024) 021 [arXiv:2401.00039] [INSPIRE].

[41] P. Dey and K. Ghosh, *Bootstrapping conformal defect operators on a line*, JHEP 10 (2024) 129 [arXiv:2404.06576] [INSPIRE].

[42] M. Serone, G. Spada and G. Villadoro, $\lambda\phi^4$ Theory I: The Symmetric Phase Beyond NNNNNNNNLO, JHEP 08 (2018) 148 [arXiv:1805.05882] [INSPIRE].

[43] M. Serone, G. Spada and G. Villadoro, $\lambda\phi^4$ theory — Part II. the broken phase beyond NNNN(NNNN)LO, JHEP 05 (2019) 047 [arXiv:1901.05023] [INSPIRE].

[44] M. Kruczenski, J. Penedones and B.C. van Rees, *Snowmass White Paper: S-matrix Bootstrap*, arXiv:2203.02421 [INSPIRE].
[45] O.J. Rosten, *Fundamentals of the Exact Renormalization Group*, Phys. Rept. **511** (2012) 177
[arXiv:1003.1366] [INSPIRE].

[46] I. Balog et al., *Convergence of Nonperturbative Approximations to the Renormalization Group*, Phys. Rev. Lett. **123** (2019) 240604 [arXiv:1907.01829] [INSPIRE].

[47] J. Preskill, *Simulating quantum field theory with a quantum computer*, PoS LATTICE2018 (2018) 024 [arXiv:1811.10085] [INSPIRE].

[48] G.K. Brennen, P. Rohde, B.C. Sanders and S. Singh, *Multiscale quantum simulation of quantum field theory using wavelets*, Phys. Rev. A **92** (2015) 032315 [arXiv:1412.0750] [INSPIRE].

[49] A. Hardy et al., *Optimized Quantum Simulation Algorithms for Scalar Quantum Field Theories*, arXiv:2407.13819 [INSPIRE].

[50] M. Hogervorst, S. Rychkov and B.C. van Rees, *Truncated conformal space approach in d dimensions: A cheap alternative to lattice field theory?*, Phys. Rev. D **91** (2015) 025005 [arXiv:1409.1581] [INSPIRE].

[51] S. Rychkov and L.G. Vitale, *Hamiltonian truncation study of the $\phi^4$ theory in two dimensions*, Phys. Rev. D **91** (2015) 085011 [arXiv:1412.3460] [INSPIRE].

[52] J. Elias-Miró, S. Rychkov and L.G. Vitale, *High-Precision Calculations in Strongly Coupled Quantum Field Theory with Next-to-Leading-Order Renormalized Hamiltonian Truncation*, JHEP **10** (2017) 213 [arXiv:1706.06121] [INSPIRE].

[53] J. Elias-Miró, S. Rychkov and L.G. Vitale, *NLO Renormalization in the Hamiltonian Truncation*, Phys. Rev. D **96** (2017) 065024 [arXiv:1706.09929] [INSPIRE].

[54] A. Tilloy, *Relativistic continuous matrix product states for quantum fields without cutoff*, Phys. Rev. D **104** (2021) 096007 [arXiv:2102.07741] [INSPIRE].

[55] A. Tilloy, *Variational method in relativistic quantum field theory without cutoff*, Phys. Rev. D **104** (2021) L091904 [arXiv:2102.07733] [INSPIRE].

[56] A. Tilloy, *A study of the quantum Singh-Gordon model with relativistic continuous matrix product states*, arXiv:2209.05341 [INSPIRE].

[57] J. Haegeman et al., *Time-Dependent Variational Principle for Quantum Lattices*, Phys. Rev. Lett. **107** (2011) 070601 [arXiv:1103.0936] [INSPIRE].

[58] L. Vanderstraeten, J. Haegeman and F. Verstraete, *Tangent-space methods for uniform matrix product states*, SciPost Phys. Lect. Notes **7** (2019) 1 [arXiv:1810.07006] [INSPIRE].

[59] F. Pollmann, S. Mukerjee, A.M. Turner and J.E. Moore, *Theory of Finite-Entanglement Scaling at One-Dimensional Quantum Critical Points*, Phys. Rev. Lett. **102** (2009) 255701 [INSPIRE].

[60] B. Vanhecke, F. Verstraete and K. Van Acoleyen, *Entanglement scaling for $\lambda\phi^2$, Phys. Rev. D **106** (2022) L071501 [arXiv:2104.10564] [INSPIRE].

[61] J. Glimm, *Boson fields with nonlinear selfinteraction in two dimensions*, Commun. Math. Phys. **8** (1968) 12.

[62] P. Federbush, *Partially Alternate Derivation of a Result of Nelson*, J. Math. Phys. **10** (1969) 50.

[63] D. Schaich and W. Loinaz, *An improved lattice measurement of the critical coupling in $\phi_2^4$ theory*, Phys. Rev. D **79** (2009) 056008 [arXiv:0902.0045] [INSPIRE].

[64] S. Bronzin, B. De Palma and M. Guagnelli, *New Monte Carlo determination of the critical coupling in $\phi_2^4$ theory*, Phys. Rev. D **99** (2019) 034508 [arXiv:1807.03381] [INSPIRE].
[65] P. Bosetti, B. De Palma and M. Guagnelli, Monte Carlo determination of the critical coupling in $\phi_2^4$ theory, Phys. Rev. D 92 (2015) 034509 [arXiv:1506.08587] [INSPIRE].

[66] A. Milsted, J. Haegeman and T.J. Osborne, Matrix product states and variational methods applied to critical quantum field theory, Phys. Rev. D 88 (2013) 085030 [arXiv:1302.5582] [INSPIRE].

[67] B. Vanhecke et al., Scaling Hypothesis for Matrix Product States, Phys. Rev. Lett. 123 (2019) 250604 [arXiv:1907.08603] [INSPIRE].

[68] D. Kadoh et al., Tensor network analysis of critical coupling in two dimensional $\phi^4$ theory, JHEP 05 (2019) 184 [arXiv:1811.12376] [INSPIRE].

[69] C. Delcamp and A. Tilloy, Computing the renormalization group flow of two-dimensional $\phi^4$ theory with tensor networks, Phys. Rev. Res. 2 (2020) 033278 [arXiv:2003.12993] [INSPIRE].

[70] J. Haegeman, J.I. Cirac, T.J. Osborne and F. Verstraete, Calculus of continuous matrix product states, Phys. Rev. B 88 (2013) 085118 [arXiv:1211.3935] [INSPIRE].

[71] J.D. Bjorken and S.D. Drell, Relativistic Quantum Mechanics, McGraw-Hill, New York, (1965) [INSPIRE].

[72] M.E. Peskin and D.V. Schroeder, An introduction to quantum field theory, Addison-Wesley, Reading, U.S.A. (1995) [DOI:10.1201/9780429503559] [INSPIRE].

[73] G.P. Lepage, Adaptive multidimensional integration: VEGAS enhanced, J. Comput. Phys. 439 (2021) 110386 [arXiv:2009.05112] [INSPIRE].

[74] M. Oshikawa and I. Affleck, Boundary conformal field theory approach to the critical two-dimensional Ising model with a defect line, Nucl. Phys. B 495 (1997) 533 [cond-mat/9612187] [INSPIRE].