text:
# Study of the 2d Ising Model with Mixed Perturbation

## Introduction

In the past the 2$`d`$ Ising model has been the subject of both analytical and numerical study. After the original solution due to Onsager  much work has been done to improve our knowledge of the model.

It has been shown that at criticality the model can be described by a Minimal Unitary Conformal Field Theory . Both magnetic and thermal perturbations of this model have been investigated and it has been found that these are the only two integrable perturbations of the model. The integrability of the model with a magnetic perturbation has been exploited by Zamolodchikov in order to obtain the exact $`S`$-matrix and mass spectrum of the theory . Further, the exact knowledge of the $`S`$-matrix has been utilized to calculate the first few terms in the spectral expansion of the two point correlation functions via the exact calculation of the form factors .

On the other hand, an infrared safe short distance expansion (IRS) for the correlators has been proposed in  (see  for a comparison of the expansions to Monte Carlo simulations). An interesting feature of the latter approach is that, in principle, it does not require the integrability of the model.

Apart from the analytic results there exists an extensive numerical investigation of the model. In particular, the transfer matrix technique has been employed in a large set of investigations: the study of the magnetic perturbation of the model , even at finite temperature , the critical equation of state  and a classification of irrelevant operators which enters in the observables of the theory .

However, the knowledge of the mixed perturbation of the model is still limited. The most important contributions to the study of this regime are due to McCoy and Wu , Delfino et al.  and to Zamolodchikov and Fonseca .

This work explores this latter regime, in fact, our purpose is to investigate the mixed perturbation model in the neighborhood of the magnetic axis. We apply the method developed by  to this case. There are two key points of this method: the high precision numerical determination of the eigenvalues of the transfer matrix as a function of the couplings and the analysis of the data by means of scaling functions obtained from a CFT approach.

The original contributions of this work can be summarized in three points. Firstly, we verify that the CFT inspired scaling functions are in perfect agreement with our data, moreover we are able to estimate the leading order corrections induced by the presence of the thermal perturbation.

Secondly, we have been able to predict the existence of a term due to the stress-energy tensor in the free energy and determine its amplitude. This is an interesting result because the contributions of the stress-energy tensor have been claimed to be zero if only integrable perturbations are involved. However, this is issue is still a matter of study.

Finally we give an estimate of the corrections to the Zamolodchikov mass spectrum and we find perfect agreement (up to some acceptable error) with Delfino et al. .

This work is organized as follows. In Section 1 we review the standard definitions used in the Ising model and we define the normalizations and conventions we use. Section 2 is devoted to the explanation of the transfer matrix technique. The computation of the scaling functions is illustrated in Section 3. Our results are presented in Section 4, and finally, we draw our conclusions in Section 5. In Appendix A we summarize all of the known results on the amplitudes of the scaling functions, in Appendix B we present an example of the obtained scaling functions.

## 1 The Ising Model

In this section we review the existing results of the 2$`d`$ Ising model. In particular, we shall present the model both in its lattice form and in its continuum formulation as a field theory. We also define the observables that we shall use in the following work.

### 1.1 The Lattice Model

The Ising model in a magnetic field at an arbitrary temperature is defined by the partition function

$`Z(\beta ,h_{\mathrm{}})={\displaystyle \underset{\sigma _i=\pm 1}{}}e^{\beta _{n,m}\sigma _n\sigma _m+h_{\mathrm{}}_n\sigma _n},`$ (1)
where the spin variable $`\sigma _n`$ takes the values $`\pm 1`$; the notation $`n,m`$ represents nearest neighbor sites on the lattice; the sites are labeled by $`n=(n_o,n_1)`$ and the two sizes of the square lattice are $`L_0`$ and $`L_1`$ (they are taken to be different because our transfer matrix calculations will treat the two directions asymmetrically); the total number of sites of the lattice will be denoted as $`N=L_0L_1`$.
The coupling $`\beta `$ is the inverse of the temperature, while the magnetic perturbation is introduced by the coupling $`h_{\mathrm{}}H\beta `$, where $`H`$ is the magnetic coupling. This model undergoes a second order phase transition when $`h_{\mathrm{}}=0`$ and $`\beta `$ reaches its critical value $`\beta _c`$

$`\beta _c{\displaystyle \frac{1}{2}}\mathrm{log}(\sqrt{2}+1)=0.4406868\mathrm{}.`$ (2)
Now we define the observables which we shall consider in the following work.

* TheMagnetization per site is defined as

$`M(\beta ,h_{\mathrm{}})={\displaystyle \frac{1}{N}}{\displaystyle \frac{}{h_{\mathrm{}}}}(\mathrm{log}Z(\beta ,h_{\mathrm{}}))={\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}\sigma _i.`$ (3)
* The Magnetic Susceptibility is defined as

$`\chi (\beta ,h_{\mathrm{}})={\displaystyle \frac{M(\beta ,h_{\mathrm{}})}{h_{\mathrm{}}}}.`$ (4)
* The Free Energy is defined as

$`f(\beta ,h_{\mathrm{}})={\displaystyle \frac{1}{N}}\mathrm{log}Z(\beta ,h_{\mathrm{}}).`$ (5)
The free energy is composed of a “bulk” term $`f_b(\beta ,h_{\mathrm{}})`$, which is an analytic function of $`h_{\mathrm{}}`$ and $`\beta `$, and a “singular” term encoding the relevant information about the theory in the neighborhood of the critical point. From the exact solution of the lattice model at $`h_{\mathrm{}}=0`$ and $`\beta =\beta _c`$ we can compute the value of $`f_b(0,0)`$

$`f_b={\displaystyle \frac{2G}{\pi }}+{\displaystyle \frac{1}{2}}\mathrm{log}2=0.9296953982\mathrm{},`$ (6)
where $`G`$ is the Catalan constant.
* The Internal Energy density is defined as

$`\widehat{E}(\beta ,h_{\mathrm{}})={\displaystyle \frac{1}{2N}}{\displaystyle \underset{n,m}{}}\sigma _n\sigma _m.`$ (7)
As for the free energy, we have both a bulk and a singular part. If we define $`E_b(\beta ,h_{\mathrm{}})`$ as the bulk contribution, we can obtain the value of $`ϵ_b=E_b(0,0)=\frac{1}{\sqrt{2}}`$ (by using Kramers-Wannier duality or CFT techniques) and we can also define (for future convenience)

$`E(\beta ,h_{\mathrm{}})=\widehat{E}(\beta ,h_{\mathrm{}})ϵ_b.`$ (8)
* Time slice correlation functions.
We can define the zero momentum projections of the two-point correlation functions $`\sigma (r)\sigma (0)`$ and $`ϵ(r)ϵ(0)`$ (they are also called time slice correlators). The magnetization of a row of the lattice (time slice) is given by

$`S_{n_0}={\displaystyle \frac{1}{L_1}}{\displaystyle \underset{n_1}{}}\sigma _{(n_0,n_1)}`$ (9)
hence, the correlation function between time slices is given by

$`G_{\sigma \sigma }^0(n_\tau )={\displaystyle \underset{n_0}{}}\left[S_{n_0}S_{n_0+n_\tau }S_{n_0}^2\right]`$ (10)
where label $`0`$ indicates that it is the zero momentum projection of the original correlator.

### 1.2 The Ising Field Theory

The 2$`d`$ Ising model near the phase transition can be described via a minimal unitary Conformal Field Theory (with central charge $`c=1/2`$) perturbed with the energy and magnetization densities $`ϵ(x)`$ and $`\sigma (x)`$

$`𝒜=𝒜_{CFT}+h{\displaystyle \text{d}^2x\sigma (x)}+\tau {\displaystyle \text{d}^2xϵ(x)}.`$ (11)
Where $`𝒜_{CFT}`$ is the action of the model at the critical temperature without external magnetic field.
The coupling constant $`h`$ represents the magnetic perturbation of the model and, as stated before, it is the continuum version of the previously defined coupling $`h_{\mathrm{}}`$. The other coupling constant $`\tau `$ represents the thermal perturbation and, near criticality, it is proportional to the reduced temperature $`t`$

$`t={\displaystyle \frac{\beta _c\beta }{\beta _c}}`$ (12)
which we use in the rest of the paper.
In the following Section we give a brief report of the known results on the field theory description of the model.

#### 1.2.1 The Critical Theory $`(t=0,h_{\mathrm{}}=0)`$

The Ising model is the lowest model of the so-called “Minimal Unitary” series of conformal theories whose central charge is given by

$`c_p=1{\displaystyle \frac{6}{p(p+1)}},p=\mathrm{\hspace{0.33em}3},\mathrm{}.`$ (13)
The peculiarity of these models is that they possess a finite set of primary fields; as a consequence the whole space of local operators of the theory can be built by applying the generators of the Virasoro algebra to the primary fields. The operators obtained in this way are called secondary fields or descendants. Following this route, one is led to organize the operator content of the theory in conformal families, i.e. the sets of descendants of each primary field.
The operator spectrum of the Ising model consists of three primary operators

* Identity $``$ $`\mathrm{\Delta }_I=0`$
* Magnetization $``$ $`\mathrm{\Delta }_\sigma =1/16`$
* Energy $``$ $`\mathrm{\Delta }_ϵ=1/2`$

where $`X_i=2\mathrm{\Delta }_i`$ is the scaling dimension of each operator.

We assume for the field $`\sigma (x)`$ and $`ϵ(x)`$ the usual CFT normalization

$`\sigma (x)\sigma (0)=|x|^{1/4},ϵ(x)ϵ(0)=|x|^2`$ (14)
so they have scaling dimensions $`1/16`$ and $`1/2`$ respectively; the parameters $`h`$ and $`t`$ have dimensions $`15/8`$ and $`1`$.
Hence there are three conformal families that descend from these operators and the general expression for the descendants $`𝒪_{[\varphi ]}`$ is

$`𝒪_{[\varphi ]}=L_{k_s\mathrm{}}L_{k_1}\overline{L}_{m_p\mathrm{}}\overline{L}_{m_1}\varphi `$ (15)
where

$`{\displaystyle \underset{i=1}{\overset{s}{}}}k_i=n;{\displaystyle \underset{i=1}{\overset{p}{}}}m_i=\overline{n},`$ (16)
the operator $`\varphi `$ is one of the primary fields of the theory.
The scaling dimension and the conformal spin of the operators are given by

$$\begin{array}{cc}\hfill X_𝒪& =\mathrm{\Delta }_𝒪+\overline{\mathrm{\Delta }}_𝒪=n+\overline{n}+\mathrm{\Delta }_\varphi +\overline{\mathrm{\Delta }}_\varphi \hfill \\ \hfill s_𝒪& =n\overline{n}.\hfill \end{array}$$
(17)
Among secondary fields, the quasi-primary fields play a special role. A descendant field $`Q`$ is called quasi-primary when

* $`L_1|Q=0`$
* $`|Q`$ is not a null vector.

It will be shown in Section 3 that quasi-primary fields can be used as building blocks for the construction of an effective Hamiltonian for the lattice model near criticality.

#### 1.2.2 The Ising Field Theory with Magnetic Perturbation $`(t=0,h_{\mathrm{}}0)`$

If $`\beta `$ is fixed to its critical value $`\beta _c=\frac{1}{2}\mathrm{log}(1+\sqrt{2})`$ and the magnetic field is switched on, the field theory is still integrable, i.e. it possesses an infinite number of integrals of motion. This implies the exact knowledge of the $`S`$-matrix and the mass spectrum. The latter result, due to Zamolodchikov , suggests that the model could be described by a scattering theory with a spectrum of eight self-conjugated particles with masses

$$\begin{array}{cc}\hfill m_2& =2m_1\mathrm{cos}\frac{\pi }{5}=(1.6180339887..)m_1,\hfill \\ \hfill m_3& =2m_1\mathrm{cos}\frac{\pi }{30}=(1.9890437907..)m_1,\hfill \\ \hfill m_4& =2m_2\mathrm{cos}\frac{7\pi }{30}=(2.4048671724..)m_1,\hfill \\ \hfill m_5& =2m_2\mathrm{cos}\frac{2\pi }{15}=(2.9562952015..)m_1,\hfill \\ \hfill m_6& =2m_2\mathrm{cos}\frac{\pi }{30}=(3.2183404585..)m_1,\hfill \\ \hfill m_7& =4m_2\mathrm{cos}\frac{\pi }{5}\mathrm{cos}\frac{7\pi }{30}=(3.8911568233..)m_1,\hfill \\ \hfill m_8& =4m_2\mathrm{cos}\frac{\pi }{5}\mathrm{cos}\frac{2\pi }{15}=(4.7833861168..)m_1\hfill \end{array}$$
(18)
where $`m_1m_1(h)`$ is the fundamental mass of the theory, and is given by

$`m_1(h)=𝒞h^{\frac{8}{15}}.`$ (19)
The numerical value of $`𝒞`$ was computed by Fateev

$`𝒞={\displaystyle \frac{4\mathrm{sin}\frac{\pi }{5}\mathrm{\Gamma }(\frac{1}{5})}{\mathrm{\Gamma }(\frac{2}{3})\mathrm{\Gamma }(\frac{8}{15})}}\left({\displaystyle \frac{4\pi ^2\mathrm{\Gamma }(\frac{3}{4})\mathrm{\Gamma }^2(\frac{13}{16})}{\mathrm{\Gamma }(\frac{1}{4})\mathrm{\Gamma }^2(\frac{3}{16})}}\right)^{\frac{4}{5}}=4.40490858\mathrm{}.`$ (20)
The vacuum expectation values of energy and magnetization can be parametrized as

$`ϵ=A_ϵh^{8/15},\sigma =A_\sigma h^{1/15}`$ (21)
where the amplitudes $`A_ϵ`$ and $`A_\sigma `$ can be computed exactly

$`A_ϵ=2.00314\mathrm{},A_\sigma =1.27758227\mathrm{}.`$ (22)
#### 1.2.3 Non-integrable Perturbation of Ising Model $`(t0,h_{\mathrm{}}0)`$

The non-integrable perturbation of the 2$`d`$ Ising model111See, e.g.  for a detailed discussion of the model with thermal perturbation only $`(t0,h_{\mathrm{}}=0)`$. was treated in  in the framework of the Form Factors Perturbation theory.
The starting point of FFPT is to consider the mixed perturbation as the perturbation of an integrable QFT.

The action (Eq. 11) describes a one parameter family of theories labeled by the adimensional scaling variable $`\xi `$

$`\xi {\displaystyle \frac{t}{h^{\frac{8}{15}}}}.`$ (23)
In this framework it is possible to calculate (at first order) the corrections to the mass spectrum and vacuum energy of the integrable theory. Two limits have been investigated in : the magnetic perturbation of a free massive Majorana fermion ($`\xi \mathrm{}`$) and the thermal perturbation of the integrable magnetic perturbation of the critical Ising model ($`\xi 0`$).
The detailed discussion of  shows that in the former case ($`\xi \mathrm{}`$), a straightforward application of the method is not an easy task. In fact, if we are in the high temperature regime of the model, the corrections at first order vanish. A second order calculation is required but is a non-trivial computation. However, in the low temperature regime, the divergence remains at first order, and is the signal kink’s confinement.

The latter case ($`\xi 0`$), however, is more tractable and gives rise to some quantitative predictions of the mass spectrum of the model. The results are listed below

$$\begin{array}{cc}\hfill \frac{\delta _{vac}}{\delta m_1}& 0.0558\mathrm{}m_1^0,\hfill \\ \hfill \frac{\delta m_2}{\delta m_1}& 0.8616\mathrm{},\hfill \\ \hfill \frac{\delta m_3}{\delta m_1}& 1.5082\mathrm{}.\hfill \end{array}$$
(24)
The non integrable perturbation gives rise to an important quantitative difference on the masses above the threshold. In fact, the integrability of the pure magnetic perturbation prevents the creation of new particles at energies above the threshold (and from any other inelastic process). The explicit breakdown of integrability implies that particles above threshold become unstable and are expected to decay. This new feature manifests itself in the corrections to the masses, namely they develop an imaginary part. This effect has been seen explicitly in the case of the first mass above the threshold $`m_4`$. The first order correction is real and given by

$`{\displaystyle \frac{\delta m_4}{\delta m_1}}`$ $``$ $`1.1460\mathrm{}`$ (25)
while, at second order, one expects a non-zero value of $`\text{Im}m_4^2`$ (it was not computed in  because the authors’ analysis only covers the first order contribution).

### 1.3 Relation Between Lattice and Continuum Operators

It is useful to define, for future convenience, the relations between the lattice and continuum definitions of the energy and magnetization operators. Near the critical point, the simplest choices for the lattice operators are

* Spin operator

$`\sigma _l(x)=\sigma _x`$ (26)
where the index $`l`$ indicates that it is a lattice discretization of the continuum operator. The magnetization per site is defined as

$`\sigma _l={\displaystyle \frac{1}{N}}{\displaystyle \underset{x}{}}\sigma _l(x).`$ (27)
* Energy operator

$`ϵ_l={\displaystyle \frac{1}{4}}\sigma _x\left[{\displaystyle \underset{y\{n.n.x\}}{}}\sigma _y\right]ϵ_b`$ (28)
where the sum runs over the nearest neighbor sites of $`x`$ and $`ϵ_b`$ is the bulk term. The energy per site is defined as

$`ϵ_l={\displaystyle \frac{1}{N}}{\displaystyle \underset{x}{}}ϵ_l(x).`$ (29)

Off-critical corrections to scaling of these operators will be discussed later.

### 1.4 Converting between Lattice and Continuum Units

In order to fix the conversions between lattice and continuum units, we shall follow the careful discussion of . One can write the lattice versions of $`\sigma `$ and $`ϵ`$ as

$$\begin{array}{cc}\hfill \sigma _l& =f_0^\sigma (t,h)\sigma +f_i(t,h)\varphi _i\hfill \\ \hfill ϵ_l& =g_0^ϵ(t,h)ϵ+g_i(t,h)\varphi _i\hfill \end{array}$$
(30)
where $`f_0^\sigma `$, $`f_i`$, $`g_0^ϵ`$, $`g_i`$ are suitable functions of $`t`$ and $`h`$ (which also depend on the parity properties of the operator), while the operators $`\varphi _i`$ are all other fields (relevant and irrelevant) of the theory respecting the symmetries of the lattice. We also have

$`h_{\mathrm{}}=b_0(t,h)h`$ (31)
which is the relation between the lattice coupling constant $`h_{\mathrm{}}`$ and the continuum magnetic field $`h`$; $`b_0(t,h)`$ is an even function of $`h`$.

At first order in $`t`$ and $`h`$, i.e. near criticality when $`t0`$ and $`h0`$, we have

$`\sigma _l=R_\sigma \sigma ,ϵ_l=R_ϵϵ,h_{\mathrm{}}=R_hh`$ (32)
where the constants $`R_\sigma `$, $`R_ϵ`$, $`R_h`$ are defined as

$$\begin{array}{cc}\hfill R_\sigma & =\underset{t,h0}{lim}f_0^\sigma (t,h)\hfill \\ \hfill R_ϵ& =\underset{t,h0}{lim}f_0^ϵ(t,h)\hfill \\ \hfill R_h& =\underset{t,h0}{lim}b_0(t,h).\hfill \end{array}$$
(33)
The previous normalizations were fixed in  by comparison with the explicit expression of the spin-spin and energy-energy critical correlation functions on the lattice. The numerical results are

$$\begin{array}{cc}\hfill R_\sigma & =0.83868\mathrm{}\hfill \\ \hfill R_ϵ& =\frac{1}{\pi }=0.31831\mathrm{}\hfill \\ \hfill R_h& =R_\sigma ^1=1.1923\mathrm{}.\hfill \end{array}$$
(34)
## 2 Transfer Matrix Technique

We face the problem of computing mass spectrum and observables by numerical diagonalization of the transfer matrix. This technique, introduced in 1941 by Kramers and Wannier , was extensively used by Baxter to obtain analytic solutions of certain statistical mechanical models . For further details and discussions about numerical implementations of the transfer matrix see, e.g. .
The basic idea is to rewrite the Boltzmann weight by means of the transfer matrix $`T(u_i,u_j)`$

$`T(u_{n_0},u_{n_0+1})=V(u_{n_0})^{1/2}U(u_{n_0},u_{n_0+1})V(u_{n_0+1})^{1/2}`$ (35)
with

$$\begin{array}{cc}\hfill U(u_{n_0},u_{n_0+1})& =\mathrm{exp}\left(\beta \underset{n_1=1}{\overset{L_1}{}}\sigma _{(n_0,n_1)}\sigma _{(n_0+1,n_1)}\right)\hfill \\ \hfill V(u_{n_0})& =\mathrm{exp}\left(\beta \underset{n_1=1}{\overset{L_1}{}}\sigma _{(n_0,n_1)}\sigma _{(n_0,n_1+1)}+h_{\mathrm{}}\underset{n_1}{\overset{L_1}{}}\sigma _{(n_0,n_1)}\right)\hfill \end{array}$$
(36)
where $`u_{n_0}=(\sigma _{(n_0,1)},\sigma _{(n_0,2)},\mathrm{},\sigma _{(n_0,L_1)})`$ is the spin configuration at the row (time-slice) $`n_0`$. The previous position implies that the partition function becomes

$$\begin{array}{cc}\hfill Z(\beta ,h_{\mathrm{}})& =\underset{\sigma _i=\pm 1}{}e^{\beta _{n,m}\sigma _n\sigma _m+h_{\mathrm{}}_n\sigma _n}\hfill \\ & =\underset{\sigma _i=\pm 1}{}T(u_1,u_2)T(u_2,u_3)\mathrm{}T(u_{L_0},u_1)\hfill \\ & =\text{tr}T^{L_0}=\underset{i}{}\lambda _i^{L_0}\hfill \end{array}$$
(37)
where $`T`$ is a positive and symmetric $`2^{L1}\times 2^{L1}`$ matrix, whose (real) eigenvalues are the $`\lambda _i`$.

### 2.1 Observables in Transfer Matrix Formalism

The numerical computation of eigenvalues and eigenvectors of the transfer matrix enables us to compute all the observables we need, provided that we specify the values of $`h_{\mathrm{}}`$, $`t`$ and $`L_1`$. We can derive all the observables from the partition function $`Z(\beta ,h_{\mathrm{}})`$ (37). For the derivation see .

* Free energy

$`f(\beta ,h_{\mathrm{}})`$ $`=`$ $`{\displaystyle \frac{1}{L_0L_1}}\mathrm{log}Z(\beta ,h_{\mathrm{}}).`$ (38)
This expression simplifies in the limit $`L_0\mathrm{}`$, in fact the leading contribution is due to the maximum eigenvalue $`\lambda _{max}`$

$`f(\beta ,h_{\mathrm{}})`$ $``$ $`{\displaystyle \frac{1}{L_1}}\mathrm{log}\lambda _{max}.`$ (39)
* Magnetization

$`\sigma _l={\displaystyle \frac{\text{tr}ST^{L_0}}{\text{tr}T^{L_0}}}`$ (40)
where $`S(u_{n_0},u_{n_1})=\delta (u_{n_0},u_{n_1})\sigma _{(n_0,1)}`$. In the limit $`L_0\mathrm{}`$ we have

$`\sigma _l=0|S|0`$ (41)
where $`|0`$ is the eigenvector associated to the $`\lambda _{max}`$ eigenvalue.
* Energy The case of the internal energy is analogous similar to that of magnetization. In the limit $`L_0\mathrm{}`$ we have

$`ϵ_l=0|E|0`$ (42)
where the matrix $`E(u_{n_0},u_{n_1})`$ is given by

$`E(u_{n_0},u_{n_1})=\delta (u_{n_0},u_{n_1})\sigma _{(n_0,1)}\sigma _{(n_0,2)}.`$ (43)
* Correlation functions and mass spectrum In the limit $`L_0\mathrm{}`$ the time slice correlation function is defined as

$`S_0S_t={\displaystyle \underset{i}{}}\mathrm{exp}(m_i|t|)0|\stackrel{~}{S}|ii|\stackrel{~}{S}|0`$ (44)
with $`\stackrel{~}{S}=\frac{1}{L_1}\delta (u_{n_0},u_{n_1})_{n_1}\sigma _{(n_0,n_1)}`$.
The mass spectrum $`m_i`$ is given by

$`m_i=\mathrm{log}\left({\displaystyle \frac{\lambda _i}{\lambda _0}}\right)`$ (45)
where the eigenvalues are organized in decreasing order of magnitude $`\lambda _{max}\lambda _0>\lambda _1>\mathrm{}>\lambda _i>\mathrm{}`$ and $`|i`$ are the normalized eigenvectors of $`T`$.

## 3 Scaling Functions

In order to study the collected data obtained from the TM, we need to know the behavior of the measured operator as a function of the perturbation variables. The fundamental ingredient of this construction is the knowledge of the whole spectrum of operators of the theory, including the OPE between them.
Hence, the operator content of the 2$`d`$ Ising model at the critical point previously discussed, enables us to build an effective Hamiltonian for the perturbed model. As discussed in detail in , the aim of this effective Hamiltonian is not to describe the model at a scale comparable with the lattice spacing; instead it has to be considered as the Hamiltonian describing the model after a suitable number of Renormalization Group transformations, i.e. at a scale that is larger with respect to the lattice spacing.

### 3.1 Lattice Construction of the Ising Model Via CFT Operators

The main idea then is to use the whole spectrum of conformal operators, defined on the continuum, to describe the corrections to scaling (due to the lattice) in the observables of the model. In order to build this Hamiltonian explicitly, we have to take into account, in principle, the following ingredients:

* Symmetries of the model
Unlike the case of the model at the critical point, which exhibits two exact symmetries ($`_2`$ and duality), the presence of the magnetic field explicitly breaks all of them. Hence, in this case, there are no constraints coming from symmetries (it is useful to remember that this argument, in the critical case, selects only the fields belonging to the conformal family of the identity).
* Symmetries of the lattice
It is crucial to define the geometry of lattice we are using in the transfer matrix calculation. In the following we consider a square lattice. This means that the rotational symmetry of the CFT is broken down to the dihedral subgroup $`D_4`$ and also operators with spin are allowed. Hence, the residual symmetry group (rotations of integer multiples of $`\pi /2`$) implies that only operators with spin $`j=4k`$, $`k`$, can appear on the lattice (see  for a detailed discussion).
* Lattice $``$ continuum relations
Lattice operators are defined in terms of continuum operators as follows

$$\begin{array}{cc}\hfill \sigma _l& =f_0^\sigma (h_{\mathrm{}},t)\sigma +h_{\mathrm{}}f_0^ϵ(h_{\mathrm{}},t)ϵ+f_i^\sigma (h_{\mathrm{}},t)\sigma _i+h_{\mathrm{}}f_i^ϵ(h_{\mathrm{}},t)ϵ_i+h_{\mathrm{}}f_i^\mathrm{𝟏}\eta _i\hfill \\ \hfill ϵ_l& =g_0^ϵ(h_{\mathrm{}},t)ϵ+h_{\mathrm{}}g_0^\sigma (h_{\mathrm{}},t)\sigma +h_{\mathrm{}}g_i^\sigma (h_{\mathrm{}},t)\sigma _i+g_i^ϵ(h_{\mathrm{}},t)ϵ_i+g_i^\mathrm{𝟏}\eta _i\hfill \end{array}$$
(46)
where $`f`$ and $`g`$ are functions of the reduced temperature $`t`$, and they are even functions of the magnetization $`h_{\mathrm{}}`$. Furthermore, the operators appearing in the previous expressions have to be compatible with the described symmetries.

Now we are able to write the following lattice (effective) Hamiltonian

$`H_{lat}=H_{CFT}+h_{\mathrm{}}\sigma +tϵ+u_i\mathrm{\Psi }_i`$ (47)
where $`\mathrm{\Psi }_i`$ are the quasi-primary fields belonging to the whole set of conformal families of the theory with spin $`j=4k`$, $`k`$.
The least irrelevant fields which enter the expression of $`H_{lat}`$ are built starting from the following quasi-primary fields of the family of the identity

$$\begin{array}{cc}\hfill Q_2^\mathrm{𝟏}& =L_2\mathrm{𝟏}\hfill \\ \hfill Q_4^\mathrm{𝟏}& =(L_2^2\frac{3}{5}L_4)\mathrm{𝟏}\hfill \end{array}$$
(48)
where the notation $`Q_n^\eta `$ is used to denote with $`\eta `$ the conformal family and with $`n`$ the level of descent. The same can be done for all the other families (for a list of low-lying quasi-primary states, up to level 10, see ).
All the results are reported in Table 1.

### 3.2 Computation of the Scaling Functions

We are now in a position to compute the singular part of the scaling functions of thermodynamic observables, e.g. the free energy, making use of the lattice Hamiltonian we have constructed. This is achieved by starting with the partition function $`Z(h_{\mathrm{}},\beta )`$ of the lattice Hamiltonian $`H_{lat}`$. Expanding this expression, it is possible to write down a formal series expansion in the variables $`h_{\mathrm{}}`$ and $`\xi 0`$. Hence we are able to obtain the scaling functions expressions 222We report in appendix B the expression of some of the scaling functions. for the non-scaling corrections to the observables of Section 1.1. To obtain these results, some remarks are in order:

* In the (numerical) transfer matrix analysis we are interested in the limit $`\xi 0`$ (the thermal perturbation is smaller than the magnetic one), so we can consider the following expressions for the VEV of the operators

$`𝒪=A_𝒪h^{\frac{dim𝒪}{dimh}}q_𝒪(\xi ,h_{\mathrm{}})`$ (49)
where $`q_𝒪`$ is an analytic function of its arguments. In the general case this ansatz is not correct because the presence of resonances induce the appearance of logarithms in the expression of VEVs.
* In the formal expansion of the partition function $`Z(h_{\mathrm{}},\beta )`$, we find that there also appear products of the conformal fields contained in $`H_{lat}`$. The correct way to deal with these is to use the fusion rules of the conformal theory

$$\begin{array}{cc}\hfill [ϵ][ϵ]& [\mathrm{𝟏}]\hfill \\ \hfill [\sigma ][ϵ]& [\sigma ]\hfill \\ \hfill [\sigma ][\sigma ]& [\mathrm{𝟏}]+[ϵ]\hfill \end{array}$$
(50)
where the notation $`[\mathrm{}]`$ means that we are referring to the whole conformal family.

#### 3.2.1 A Peculiar Case: The Free Energy

In this Section we show how to determine the scaling functions of the model in the region of interest ($`\xi 0`$). We have computed the scaling functions by means of two different methods: the renormalization group, and CFT approach.
If we consider the free energy as an example, in the renormalization group approach, we can write it as the sum of three contributions

$`f(t,h_{\mathrm{}})=f_b(t,h_{\mathrm{}})+f_{sing}(t,h_{\mathrm{}})+f_{log}(t,h_{\mathrm{}})`$ (51)
The bulk term takes into account analytic contributions in the variables $`t`$ and $`h_{\mathrm{}}`$ due to non critical behavior

$$\begin{array}{cc}\hfill f_b(t,h_{\mathrm{}})& =f_{0,0}^b+f_{0,2}^bh_{\mathrm{}}^2+f_{0,4}^bh_{\mathrm{}}^4+(f_{1,0}^b+f_{1,2}^bh_{\mathrm{}}^2+f_{1,4}^bh_{\mathrm{}}^4)t+\hfill \\ & +(f_{2,0}^b+f_{2,2}^bh_{\mathrm{}}^2+f_{2,4}^bh_{\mathrm{}}^4)t^2+(f_{3,0}^b+f_{3,2}^bh_{\mathrm{}}^2+f_{3,4}^bh_{\mathrm{}}^4)t^3+O(h_{\mathrm{}}^6,t^4).\hfill \end{array}$$
(52)
Only even powers of $`h_{\mathrm{}}`$ appear because the free energy is even under $`_2`$ transformations.
The non-analytic contribution is given by the master equation of the RG

$`f_{sing}(t,h_{\mathrm{}})=g_h^{2/\mathrm{\Delta }}Y({\displaystyle \frac{g_t}{g_h^{1/\mathrm{\Delta }}}},\left\{g_ug_h^{|y_u|/\mathrm{\Delta }}\right\})`$ (53)
where the scaling variables $`g_h(t,h_{\mathrm{}})`$, $`g_t(t,h_{\mathrm{}})`$, $`g_u(t,h_{\mathrm{}})`$ are defined in the usual way

$$\begin{array}{cc}\hfill g_t(t,h_{\mathrm{}})& =t+b_th_{\mathrm{}}^2+c_tt^2+d_tt^3+e_tth_{\mathrm{}}^2+f_tt^4+g_th_{\mathrm{}}^4+h_tt^2h_{\mathrm{}}^2+\mathrm{}\hfill \\ \hfill g_h(t,h_{\mathrm{}})& =h_{\mathrm{}}(1+c_ht+d_ht^2+e_hh_{\mathrm{}}^2+f_ht^3+g_hth_{\mathrm{}}^2)+\mathrm{}\hfill \\ \hfill g_u(t,h_{\mathrm{}})& =u+a_ut+b_uh_{\mathrm{}}^2+c_ut^2+d_ut^3+e_uth_{\mathrm{}}^2+f_uh_{\mathrm{}}^4+l_ut^2h_{\mathrm{}}^2+\mathrm{}\hfill \end{array}$$
(54)
and some of the coefficients are known either exactly or numerically (see appendix A).
It has been known for a long time  that, in order to take into account the logarithmic divergence of the specific heat of the Onsager solution , the term $`f_{log}`$ (51) must have the following form

$`f_{log}=g_t^2(t,h_{\mathrm{}})\mathrm{log}\left(|g_t(t,h_{\mathrm{}})|^1\right)\stackrel{~}{Y}({\displaystyle \frac{g_h}{g_t^\mathrm{\Delta }}},\left\{g_ug_h^{|y_u|/\mathrm{\Delta }}\right\}).`$ (55)
The function $`\stackrel{~}{Y}(,)`$ can be considered as a constant in our analysis, whose exact value is $`\stackrel{~}{Y}(,)A_f=\frac{4\beta _c^2}{\pi }`$.
The previous expression needs some clarification in order to extend it around the magnetic axis. Following the discussion of , we point out that if we want $`f(t,h_{\mathrm{}})`$ to be an analytic function of $`t`$ at fixed $`h_{\mathrm{}}`$, a logarithmic term should appear in the expansion of $`Y(,)`$. In order to make the $`\mathrm{log}|g_t|`$ term disappear, leaving us with the correct $`\mathrm{log}|g_h|`$, we have to extract a term like

$`x^2|g_h|^{\frac{2}{\mathrm{\Delta }}}\mathrm{log}|x|^1A_f`$ (56)
from the expansion of $`Y(,)`$ when $`x0`$, where $`x`$ is the adimensional variable $`\frac{g_t}{g_h^{1/\mathrm{\Delta }}}`$. Hence the logarithmic contribution has the following form 333In principle it is possible to include higher logarithmic powers in the free energy scaling function, However the leading order of the expansion of such terms are zero with a precision of $`10^6`$ in our fits.

$`f_{log}(t,h_{\mathrm{}})=g_t^2\mathrm{log}|g_h|^1{\displaystyle \frac{A_f}{\mathrm{\Delta }}}.`$ (57)
This gives the correct divergence of the specific heat

$`C{\displaystyle \frac{^2f}{t^2}}2{\displaystyle \frac{A_f}{\mathrm{\Delta }}}\mathrm{log}|h_{\mathrm{}}|+\mathrm{}.`$ (58)
With regard to the CFT approach one remark is in order: the CFT derivation can generate only the singular term ($`f_{sing}`$),while the bulk and the logarithmic contributions must be added separately.

The comparison of the two different approaches (CFT and Renormalization Group) enables us to understand the relation between the operators and the corrections they give rise to (for a careful discussion about this crucial point see ). A relevant example is given by the stress-energy tensor: it is the only responsible for terms like $`h^{32/15}`$, $`th^{32/15}`$, $`t^2h^{32/15}`$. This will have fundamental implications in discussing the numerical data.

## 4 Analysis of Results

The analysis of the data obtained by transfer matrix technique is done as follows:

* To analyze the large amount of data, due to the large number of required terms in the expansion of the scaling functions, and in order to obtain sensible results in the fits, we fix as many terms as possible, resorting to both exact and high precision numerical results already known in the literature;
* Fixing all known parameters in the scaling function and then fitting the data enables us to check their correctness in this particular regime ($`\xi 0`$, i.e. the neighborhood of the magnetic axis);
* The high precision of our data enables us to conjecture the presence or the absence of contributions due to well identified sources, i.e. terms like $`(At+Bt^2)h^{32/15}`$, which are entirely due to the stress-energy tensor (the irrelevant fields $`T\overline{T}`$, $`T^2`$ and $`\overline{T}`$).

This Section is devoted to developing such an analysis and to discuss our results.

### 4.1 Outline of the Numerical Computations

We perform our numerical computations extracting the four larger eigenvalues of the transfer matrix at given $`t`$, $`h_{\mathrm{}}`$, and $`L`$, the size of the transverse direction of the lattice.
The diagonalization method utilizes an iterative algorithm, due to , that evaluates the highest-lying eigenvalues of the transfer matrix with arbitrary high precision. The choice of an iterative algorithm lies in the reason that the typical size of the reduced transfer matrix is of order $`10^5`$ and the exact diagonalization, even if partial, is an unobtainable computationally.
Our programs we run on 10 PCs equipped with Pentium III processors and 256 Mb of RAM for a total CPU time of about 2 months.
We performed our numerical computations for each available choice of the three parameters $`t`$, $`h_{\mathrm{}}`$ and $`L`$ that we collected in table 2 for a total of 2600 different runs.

### 4.2 Infinite Volume Extrapolations

After obtaining the values of the thermodynamical quantities at fixed $`L`$, we would like to extrapolate them, taking the thermodynamic limit. In order to do this, we follow the method of , which we outline briefly for sake of completeness.
The basic idea is that the behavior of any thermodynamic quantity in a massive QFT will show an exponential decay as a function of $`L`$. The task we would like to achieve is to find the asymptotic value of the observables; in order to reach this result, we must be able to subtract all the possible exponential behaviors from our data. To achieve we iteratively subtract the exponential behavior from our data, removing both leading and subleading corrections. We define

$$\begin{array}{cc}\hfill b^i(L_12)& =c\mathrm{exp}(x(L_12))+b^{i+1}(L_1)\hfill \\ \hfill b^i(L_11)& =c\mathrm{exp}(x(L_11))+b^{i+1}(L_1)\hfill \\ \hfill b^i(L_1)& =c\mathrm{exp}(x(L_1))+b^{i+1}(L_1)\hfill \end{array}$$
(59)
where $`b^0(L_1)`$ is the transfer matrix quantity. A step of iteration is defined by solving the previous system with respect to $`b^{i+1}(L_1)`$, $`c`$ and $`x`$. The iteration chain stops when the predicted value becomes numerically unstable. The taken result is the last stable prediction and its error is evaluated from the variation with respect to previous step.

### 4.3 The Fitting Procedure

In order to analyze our data we fit them with the scaling function previously obtained for a generic operator $`𝒪`$.
The fit is performed in three steps:

1. We start fitting the data at fixed $`h_{\mathrm{}}`$ with a polynomial in the reduced temperature

$`𝒪(t)q|_h_{\mathrm{}}𝒪(0)|_h_{\mathrm{}}=𝒜_1^𝒪(h_{\mathrm{}})t+𝒜_2^𝒪(h_{\mathrm{}})t^2+𝒜_3^𝒪(h_{\mathrm{}})t^3+𝒜_4^𝒪(h_{\mathrm{}})t^4+\mathrm{}`$ (60)
where we found convenient to subtract the value of the observable at the critical temperature; with this subtraction we are a position to discard all the terms in the expansion that depend on $`h_{\mathrm{}}`$ only. The results of the fit are collected in a table displaying the values of $`𝒜_1^𝒪(h_{\mathrm{}})`$ and $`𝒜_2^𝒪(h_{\mathrm{}})`$ as a functions of $`h_{\mathrm{}}`$.
2. By means of the expansion for the scaling function, we can fit the functions $`𝒜_1^𝒪(h_{\mathrm{}})`$ and $`𝒜_2^𝒪(h_{\mathrm{}})`$ against the magnetic field $`h_{\mathrm{}}`$. In this way we were able to obtain the amplitudes of the corrections to the observable $`𝒪`$ at first and second order in the reduced temperature $`t`$.
3. We use all predicted values of the parameters of the fit to calculate the $`\chi ^2`$ with a unique fitting function for all the data (the coefficients $`𝒜_3^𝒪(h_{\mathrm{}})`$ and $`𝒜_4^𝒪(h_{\mathrm{}})`$ of step 1 are only useful to the determination of the $`\chi ^2`$ of the entire set of data).

The list of requirements a fit must fulfill to be accepted is the following:

* The reduced $`\chi ^2`$ of the fit had to be of order unity, i.e. we required it to have a confidence level larger than $`30\%`$.
* In order to be included in the fitting function, the subleading terms must have an amplitude larger than the corresponding error.
* The number of degrees of freedom must be larger than 5.

The previous requirements are very hard to fulfill at the same time; this way we also take into account the systematic errors due to extrapolations.

### 4.4 Numerical Results

#### 4.4.1 Determination of $`e_h`$

The least precise of the constants collected in the appendix A is $`e_h`$. The high precision of our data enables us to make an attempt to improve this estimate. Following the same procedure of , and fixing all the known terms in the scaling function we are able to give the following result

$`e_h=0.007298(3)`$ (61)
where the data we used as input are collected in appendix A.

#### 4.4.2 Free Energy

The analysis of the free energy can be performed in four steps.

1. we fit the the term proportional to $`t`$ with the following function

$$\begin{array}{cc}& 0.6232260.511645h_{\mathrm{}}^{\frac{8}{15}}+0.329993h_{\mathrm{}}^{\frac{16}{15}}0.0427280h_{\mathrm{}}^2+\hfill \\ & 0.0107535h_{\mathrm{}}^2\mathrm{log}|h_{\mathrm{}}|+h_{\mathrm{}}^{\frac{32}{15}}𝒦_0+0.0184589h_{\mathrm{}}^{\frac{38}{15}}+h_{\mathrm{}}^{\frac{40}{15}}𝒦_1+\hfill \\ & +0.0001605h_{\mathrm{}}^{\frac{46}{15}}+h_{\mathrm{}}^{\frac{48}{15}}𝒦_2+h_{\mathrm{}}^{\frac{52}{15}}𝒦_3+\mathrm{}.\hfill \end{array}$$
(62)
As we explain in the appendix A we set $`Y^{(2,1)}(0,0)=0`$ according to the results of  and , moreover we can found $`Y^{(1,1)}(0,0)`$ compatible with zero within the error in all our fits.
The fit shows a non-zero correction due to the stress-energy tensor, in fact we are able to estimate the value of $`𝒦_0`$

$`0.05828<𝒦_0<0.05820.`$ (63)
where the expression of $`𝒦_0`$ in terms of the scaling fields is given by

$`𝒦_0=\left(a_u+0.664773u\right)Y^{(0,1)}(0,0).`$ (64)
This is quite an interesting result, in fact it is known that in both the integrable perturbation of the critical Ising model, the amplitude of this term is compatible with zero.
2. we perform the fit of $`t^2`$ term with the following function

$$\begin{array}{cc}& 𝒦_40.131877\mathrm{log}|h_{\mathrm{}}|0.244467h_{\mathrm{}}^{\frac{8}{15}}+0.299064h_{\mathrm{}}^{\frac{16}{15}}+\hfill \\ & +𝒦_5h_{\mathrm{}}^{\frac{22}{15}}0.0318458h_{\mathrm{}}^2+0.0100527h_{\mathrm{}}^2\mathrm{log}|h_{\mathrm{}}|+h_{\mathrm{}}^{\frac{32}{15}}𝒦_6+\mathrm{}.\hfill \end{array}$$
(65)
As previously stated we put $`Y^{(2,1)}(0,0)=Y^{(1,1)}(0,0)=0`$, and we found that leading correction is

$`0.03655<𝒦_4<0.03657.`$ (66)
The expression of $`𝒦_4`$ in terms of the scaling fields is:

$`𝒦_4=0.481290866+f_{2,0}^b`$ (67)
in this way we are able to give an estimate for the bulk term

$`0.44474<f_{2,0}^b<0.44472.`$ (68)
Furthermore, we can clearly see that the term due to the stress-energy tensor $`𝒦_6`$ (which comes with the power $`h_{\mathrm{}}^{32/15}`$) is different from zero, as for the previous case.
As discussed in  it is difficult to have reliable estimates for the amplitudes of the subleading terms444As discussed in , notwithstanding systematic errors, one can give also a rough estimate of the subleading amplitudes. Following the same route we are able to find    $`0.721<𝒦_5<0.717.`$  (69)  due to large systematic deviations induced by the uncertainties in the leading corrections, nevertheless we can still assert that they are different from zero.
3. we check that within the precision of our computation, the function $`\stackrel{~}{Y}(,)`$ is a constant. In fact, if we expand it in Taylor series, we find that the first contribution to the scaling function comes in the $`t^2`$ term and has the following form

$`\stackrel{~}{Y}^{(1,0)}(0,0)\mathrm{log}|h_{\mathrm{}}|h_{\mathrm{}}^{22/15}.`$ (70)
If we perform the fit with this new contribution, we find that its amplitude is compatible with zero. Hence, for our purposes it is safe to consider $`\stackrel{~}{Y}(,)`$ as a constant.
4. we compute the $`\chi ^2`$ with the global scaling function (Eq. 60) both on $`t`$ and $`h_{\mathrm{}}`$ utilizing the value of the constants predicted in the previous steps, in order to verify the correctness of the fitting procedure.

It is important to mention that in the fitting functions we use the known values of the constants of Appendix A. It is a non-trivial test on the validity of both our results, and the well known values we used.

#### 4.4.3 Magnetization

1. we write down the correction proportional to the $`t`$ term

$$\begin{array}{cc}& \frac{0.2728775}{h_{\mathrm{}}^{\frac{7}{15}}}+0.3519930h_{\mathrm{}}^{\frac{1}{15}}0.0747026h_{\mathrm{}}+0.0215069h_{\mathrm{}}\mathrm{log}|h_{\mathrm{}}|+\hfill \\ & +h_{\mathrm{}}^{\frac{17}{15}}_0+0.0467625h_{\mathrm{}}^{\frac{23}{15}}+h_{\mathrm{}}^{\frac{5}{3}}_10.0004923h_{\mathrm{}}^{\frac{31}{15}}+h_{\mathrm{}}^{\frac{11}{5}}_2+\hfill \\ & +h_{\mathrm{}}^{\frac{37}{15}}_3+h_{\mathrm{}}^{\frac{41}{15}}_4+\mathrm{}\hfill \end{array}$$
(71)
where all the known quantities are taken into account. In particular we also impose the constraints we found in the analysis of the free energy, i.e. $`Y^{(2,1)}(0,0)=Y^{(1,1)}(0,0)=0`$, and we find perfect agreement also in this case.
The leading correction is due to the stress-energy tensor, as we expected from the analysis of the free energy, and its amplitude is given by

$`0.124<_0<0.122`$ (72)
which is compatible with $`_0\frac{32}{15}𝒦_0`$ within our errors.
2. the contribution due to $`t^2`$ is given by

$$\begin{array}{cc}& \frac{0.1318769}{h_{\mathrm{}}}\frac{0.1303825}{h_{\mathrm{}}^{\frac{7}{15}}}+0.3190016h_{\mathrm{}}^{\frac{1}{15}}+\frac{22}{15}h_{\mathrm{}}^{\frac{7}{15}}_50.0536389h_{\mathrm{}}+\hfill \\ & +0.0201054h_{\mathrm{}}\mathrm{log}|h_{\mathrm{}}|+h_{\mathrm{}}^{\frac{17}{15}}_6+\mathrm{}\hfill \end{array}$$
(73)
and is obtained following the same strategy as before. This order also shows a non-zero contribution due to stress-energy tensor, i.e. $`_6`$.
The amplitude of leading order correction is given by

$`0.723<_5<0.718.`$ (74)
We can check that it is in reasonable agreement with the subleading amplitude of the free energy $`_5𝒦_5`$ (see (69)).

#### 4.4.4 Internal Energy

1. the scaling function for the $`t`$ term is

$$\begin{array}{cc}& 𝒲_0+0.2992532\mathrm{log}|h_{\mathrm{}}|+0.5547414h_{\mathrm{}}^{\frac{8}{15}}+h_{\mathrm{}}^{\frac{16}{15}}𝒲_1+h_{\mathrm{}}^{\frac{22}{15}}𝒲_20.0918396h_{\mathrm{}}^2+\hfill \\ & +h_{\mathrm{}}^{\frac{32}{15}}𝒲_3+h_{\mathrm{}}^{\frac{38}{15}}𝒲_40.0228115h_{\mathrm{}}^2\mathrm{log}|h_{\mathrm{}}|+\mathrm{}\hfill \end{array}$$
(75)
the predicted value for the leading contribution is

$`0.1659<𝒲_0<0.1657`$ (76)
which expressed in term of scaling function is $`𝒲_0(2.1842763+4.5383706f_{2,0}^b)`$, the value of $`f_{2,0}^b`$ is consistent with the previous estimate (68). As before the fit is compatible with the constraint $`Y^{(2,1)}(0,0)=Y^{(1,1)}(0,0)=0`$.
2. the scaling function for the $`t^2`$ term is

$$\begin{array}{cc}& \frac{𝒲_5}{h_{\mathrm{}}^{\frac{8}{15}}}+𝒲_60.2797532\mathrm{log}|h_{\mathrm{}}|+h_{\mathrm{}}^{\frac{8}{15}}𝒲_7+h_{\mathrm{}}^{\frac{14}{15}}𝒲_6+h_{\mathrm{}}^{\frac{16}{15}}𝒲_8+\hfill \\ & 0.0727482h_{\mathrm{}}^{\frac{22}{15}}𝒲_5+\mathrm{}\hfill \end{array}$$
(77)
the predicted value for the leading contribution is

$`0.414<𝒲_5<0.408`$ (78)
where $`𝒲_50.5672963Y^{(3,0)}(0,0)`$, and again the fit is compatible with the constraint $`Y^{(2,1)}(0,0)=Y^{(1,1)}(0,0)=0`$.

#### 4.4.5 Susceptibility

The form of the scaling function is in total agreement with the computed values of the this observable. However no new predictions are available.

#### 4.4.6 Mass spectrum

In order to check the results of the FFPT proposed by , we performed the fit on the square of the first three masses of the theory, and computed the following ratio

$`{\displaystyle \frac{\delta m_i^2}{\delta m_1^2}}={\displaystyle \frac{m(t)_i^2m(0)_i^2}{m(t)_1^2m(0)_1^2}}`$ (79)
where $`m(t)_i^2`$ are the perturbed masses, and $`m(0)_i^2`$ the unperturbed ones.
The results of our analysis are reported below:

$$\begin{array}{cc}\hfill \frac{\delta m_2^2}{\delta m_1^2}& =1.393(17)\hfill \\ \hfill \frac{\delta m_3^2}{\delta m_1^2}& =3.16(30).\hfill \end{array}$$
(80)
These results are in perfect agreement with the theoretical prediction. In fact, since the deviation from the integrable model is small, it is correct to write $`\delta m_i^2=2m(0)_i\delta m_i`$, and finally we have

$$\begin{array}{cc}\hfill \frac{\delta m_2}{\delta m_1}& =0.86(1)\hfill \\ \hfill \frac{\delta m_3}{\delta m_1}& =1.58(15).\hfill \end{array}$$
(81)
#### 4.4.7 Magnetization Overlap

The magnetic overlap $`|F_1^\sigma |^2`$ can also be analyzed as before. We are able to find the leading order corrections of terms proportional to $`t`$

$$\begin{array}{c}\hfill \frac{_0}{h^{\frac{8}{15}}}+h^{\frac{2}{15}}_1+h^{\frac{8}{15}}_2+h^{\frac{14}{15}}_3+h^{\frac{16}{15}}_4+h^{\frac{22}{15}}_5+h^{\frac{8}{5}}_6+h^{\frac{28}{15}}_7+h^2_8+\mathrm{}\end{array}$$
(82)
and $`t^2`$

$$\begin{array}{c}\hfill \frac{_9}{h^{\frac{16}{15}}}+\frac{_{10}}{h^{\frac{8}{15}}}+\frac{_{11}}{h^{\frac{2}{15}}}+h^{\frac{2}{5}}_{12}+h^{\frac{8}{15}}_{13}+h^{\frac{14}{15}}_{14}+\mathrm{}.\end{array}$$
(83)
The numerical values are

$$\begin{array}{c}\hfill 0.628<_0<0.631\\ \hfill 0.661<_9<0.664\end{array}$$
(84)
It would be interesting to calculate the same corrections on theoretical grounds (at least for the $`t`$ correction) in order to make a comparison.

## 5 Discussion and Conclusions

In this paper we studied the effect of a mixed relevant perturbation on the Ising model using the Transfer Matrix technique. We consider the neighborhood of the magnetic axis in the limit $`\xi 1`$, where $`\xi `$ is the adimensional parameter $`t/h^{8/15}`$.

We have concentrated our efforts on the following areas

* We calculated the scaling functions (appendix B) applying a CFT approach to the problem.
* We used all known predictions about the behavior of scaling functions (see appendix A), and verified that they all agree with our data.
* Being able to identify the contribution of secondary fields to the scaling function, we predicted that there is a non zero contribution due to the stress-energy tensor, and we evaluate it. In our opinion this a quite interesting result, because it is known that the contribution of these particular secondary fields is zero if we study the model with only one relevant perturbation.
* We obtained estimates of several amplitudes never predicted by any other analytical method before.
* We calculated the correction to the Zamolodchikov mass spectrum of the Ising model with a magnetic field, and we found perfect agreement with Delfino et al. .

There are two possible developments of this work: We can use our data to improve the knowledge of the equation of state of the Ising model, in order to map all the possible regime of perturbations. We can extend the analysis of  to study the effect of mixed perturbation on the finite temperature results.

## Acknowledgments

We thank M.Caselle and H.Hasenbusch for useful suggestions, and several useful discussions.

## Appendix A Known Numbers

It is known that it is possible to write the free energy of the model in terms of nonlinear scaling fields .

The scaling fields are analytic functions of $`t`$ and $`h_{\mathrm{}}`$ respecting the $`_2`$ parity of $`h_{\mathrm{}}`$. Their Taylor expansions are expected to be

$$\begin{array}{cc}\hfill g_t(t,h_{\mathrm{}})& =t+b_th_{\mathrm{}}^2+c_tt^2+d_tt^3+e_tth_{\mathrm{}}^2+f_tt^4+g_th_{\mathrm{}}^4+h_tt^2h_{\mathrm{}}^2+\mathrm{}\hfill \\ \hfill g_h(t,h_{\mathrm{}})& =h_{\mathrm{}}(1+c_ht+d_ht^2+e_hh_{\mathrm{}}^2+f_ht^3+g_hth_{\mathrm{}}^2)+\hfill \\ \hfill g_u(t,h_{\mathrm{}})& =u+a_ut+b_uh_{\mathrm{}}^2+c_ut^2+d_ut^3+e_uth_{\mathrm{}}^2+f_uh_{\mathrm{}}^4+l_ut^2h_{\mathrm{}}^2+\mathrm{}\hfill \end{array}$$
(85)
Here we report all the analytically known coefficients

$$c_h=\frac{\beta _c}{\sqrt{2}},d_h=\frac{23\beta _c^2}{16},f_h=\frac{191\beta _c^3}{48\sqrt{2}},$$

$$c_t=\frac{\beta _c}{\sqrt{2}},d_t=\frac{7\beta _c^2}{6},f_t=\frac{17\beta _c^3}{6\sqrt{2}},$$
(86)
$$e_t=b_t\beta _c\sqrt{2},b_t=\frac{E_0\pi }{16\beta _c^2}$$

where

$`E_0=0.0403255003\mathrm{}`$ (87)
From the analysis of the model at the critical temperature and $`h_{\mathrm{}}0`$, we obtain a new estimate for $`e_h`$

$`e_h=0.007298(3).`$ (88)
For the free energy, we use also the coefficient reported in the following

$$\begin{array}{cc}\hfill Y(0,0)& =0.99279949\mathrm{}\hfill \\ \hfill Y^{(1,0)}(0,0)& =A_{f,2}^l/b_t=0.511645336\mathrm{}\hfill \\ \hfill A_{f,2}^l& =\frac{A_E^l\pi }{Y(0,0)8\beta _c}=0.0208602\mathrm{}\hfill \\ \hfill f_b& =\frac{2}{\pi }G+\frac{1}{2}\mathrm{log}2\hfill \\ \hfill A_{f,b}^l& =0.0524442\mathrm{}\hfill \end{array}$$
(89)
where $`G`$ is the Catalan constant and the $`A`$’s are defined in .
A high precision study of the thermal perturbation have been performed by Orrick et al. in , from this work we are able to extract a set of parameters for our scaling functions.
In particular we are able to observe that the known coefficients for both thermal and magnetic perturbations are exactly consistent, and we can extend the knowledge about the thermal coefficients to gain predictions about the unknown magnetic terms.
To achieve these results we compare the scaling function of the susceptibility along the thermal axis with the high precision estimates of  and , in this way we obtained

$$\begin{array}{c}\hfill Y^{(2,0)}(0,0)=0.9625817322\mathrm{}f_{0,2}^bA_{f,b}^lY(0,0)=0.05206662255\mathrm{}f_{1,2}^b=0.00348278\mathrm{}f_{2,2}^b=0.000528775\mathrm{}.\end{array}$$
(90)
These results are obtained for the thermal perturbation in the regime $`t0,h_{\mathrm{}}=0`$, we observe that these results are valid also in our regime of interest, because the analytic continuation of Section 3.2.1 do not affect similar kinds of terms, in total agreement with our fits.
Furthermore we obtain also the relations:

$$\begin{array}{c}\hfill Y^{(2,1)}(0,0)(a_u+0.701128u)=0\\ \hfill Y^{(2,1)}(0,0)u=0\end{array}$$
(91)
In order to fulfill the above requirements, we set $`Y^{(2,1)}(0,0)=0`$ because the other choice $`u=0`$, $`a_u=0`$ was not consistent with our fits.

Finally, from Onsager’s exact solution  we shall extract

$`f_{1,0}^b=0.623226\mathrm{}.`$ (92)
## Appendix B Scaling Functions

Here we report some of the scaling functions we obtained. We remark that the functions have been evaluated to higher order of the expansion reported here ($`O(h_{\mathrm{}}^5)O(t^5)`$).
Free Energy:

$`f(t)|_h_{\mathrm{}}f(0)|_h_{\mathrm{}}=𝒜_1^f(h_{\mathrm{}})t+𝒜_2^f(h_{\mathrm{}})t^2+\mathrm{}`$ (93)
$$\begin{array}{cc}\hfill 𝒜_1^f(h_{\mathrm{}})& =f_{1,0}^b+Y^{(1,0)}(0,0)h_{\mathrm{}}^{\frac{8}{15}}+\frac{16}{15}c_hY(0,0)h_{\mathrm{}}^{\frac{16}{15}}+uY^{(1,1)}(0,0)h_{\mathrm{}}^{\frac{8}{5}}+\hfill \\ & +\left(\frac{16}{15}A_f\mathrm{log}|h_{\mathrm{}}|b_t+f_{1,2}^b+b_tY^{(2,1)}(0,0)\right)h_{\mathrm{}}^2+\mathrm{}\hfill \end{array}$$
(94)
$$\begin{array}{cc}\hfill 𝒜_2^f(h_{\mathrm{}})& =(\frac{8}{15}A_f\mathrm{log}|h_{\mathrm{}}|+f_{0,2}^b+\frac{1}{2}Y^{(2,0)}(0,0))+(\frac{8}{15}c_hY^{(1,0)}(0,0)+\hfill \\ & +c_tY^{(1,0)}(0,0))h_{\mathrm{}}^{\frac{8}{15}}+(\frac{8}{225}c_{h}^{}{}_{}{}^{2}Y(0,0)+\frac{16}{15}d_hY(0,0)+\frac{1}{2}uY^{(2,1)}(0,0))h_{\mathrm{}}^{\frac{16}{15}}+\hfill \\ & +\frac{1}{2}b_tY^{(3,0)}(0,0)h_{\mathrm{}}^{\frac{22}{15}}+\left(a_uY^{(1,1)}(0,0)+\frac{8}{5}uc_hY^{(1,1)}(0,0)+uc_tY^{(1,1)}(0,0)\right)h_{\mathrm{}}^{\frac{8}{5}}+\hfill \\ & +(\frac{16}{15}A_fb_tc_h\frac{16}{15}A_f\mathrm{log}|h_{\mathrm{}}|b_tc_t\frac{8}{15}A_fe_h\frac{16}{15}A_f\mathrm{log}|h_{\mathrm{}}|e_t+f_{2,2}^b+\hfill \\ & +b_tc_tY^{(2,0)}(0,0)+e_tY^{(2,0)}(0,0))h_{\mathrm{}}^2+\mathrm{}\hfill \end{array}$$
(95)
Internal Energy:

$`E(t)|_h_{\mathrm{}}E(0)|_h_{\mathrm{}}=𝒜_1^E(h_{\mathrm{}})t+𝒜_2^E(h_{\mathrm{}})t^2+\mathrm{}`$ (96)
$$\begin{array}{cc}\hfill 𝒜_1^E(h_{\mathrm{}})& =\left(\frac{16}{15}A_f\mathrm{log}|h_{\mathrm{}}|T_c+2f_{2,0}^bT_c+T_cY^{(2,0)}(0,0)\right)+\hfill \\ & +\left(\frac{16}{15}c_hT_cY^{(1,0)}(0,0)+2c_tT_cY^{(1,0)}(0,0)\right)h_{\mathrm{}}^{\frac{8}{15}}+\hfill \\ & +\left(\frac{16}{225}c_{h}^{}{}_{}{}^{2}T_cY(0,0)+\frac{32}{15}d_hT_cY(0,0)+uT_cY^{(2,1)}(0,0)\right)h_{\mathrm{}}^{\frac{16}{15}}+\hfill \\ & +b_tT_cY^{(3,0)}(0,0)h_{\mathrm{}}^{\frac{22}{15}}+(2a_uT_cY^{(1,1)}(0,0)+\hfill \\ & +\frac{16}{5}uc_hT_cY^{(1,1)}(0,0)+2uc_tT_cY^{(1,1)}(0,0))h_{\mathrm{}}^{\frac{8}{5}}+\mathrm{}\hfill \end{array}$$
(97)
$$\begin{array}{cc}\hfill 𝒜_2^E(h_{\mathrm{}})& =\frac{T_c}{2h_{\mathrm{}}^{\frac{8}{15}}}Y^{(3,0)}(0,0)+(\frac{8}{5}A_fc_hT_c\frac{16}{5}A_f\mathrm{log}|h_{\mathrm{}}|c_tT_c+3f_{3,0}^bT_c+\hfill \\ & +\mathrm{\hspace{0.17em}3}c_tT_cY^{(2,0)}(0,0))+(\frac{28}{75}c_{h}^{}{}_{}{}^{2}T_cY^{(1,0)}(0,0)+\frac{8}{5}c_hc_tT_cY^{(1,0)}(0,0)+\hfill \\ & +\frac{8}{5}d_hT_cY^{(1,0)}(0,0)+3d_tT_cY^{(1,0)}(0,0)+\frac{1}{2}uT_cY^{(3,1)}(0,0))h_{\mathrm{}}^{\frac{8}{15}}+\hfill \\ & +\frac{1}{2}b_tT_cY^{(4,0)}(0,0)h_{\mathrm{}}^{\frac{14}{15}}+(\frac{112}{3375}c_{h}^{}{}_{}{}^{3}T_cY(0,0)+\frac{16}{75}c_hd_hT_cY(0,0)+\hfill \\ & +\frac{16}{5}f_hT_cY(0,0)+\frac{3}{2}a_uT_cY^{(2,1)}(0,0)+\frac{8}{5}uc_hT_cY^{(2,1)}(0,0)+\hfill \\ & +3uc_tT_cY^{(2,1)}(0,0))h_{\mathrm{}}^{\frac{16}{15}}+(\frac{4}{5}b_tc_hT_cY^{(3,0)}(0,0)+3b_tc_tT_cY^{(3,0)}(0,0)+\hfill \\ & \frac{4}{15}e_hT_cY^{(3,0)}(0,0)+\frac{3}{2}e_tT_cY^{(3,0)}(0,0))h_{\mathrm{}}^{\frac{22}{15}}+(\frac{24}{5}a_uc_hT_cY^{(1,1)}(0,0)+\hfill \\ & +\frac{36}{25}uc_{h}^{}{}_{}{}^{2}T_cY^{(1,1)}(0,0)+3a_uc_tT_cY^{(1,1)}(0,0)+\frac{24}{5}uc_hc_tT_cY^{(1,1)}(0,0)+\hfill \\ & +3c_uT_cY^{(1,1)}(0,0)+\frac{24}{5}ud_hT_cY^{(1,1)}(0,0)+3ud_tT_cY^{(1,1)}(0,0)+\hfill \\ & +\frac{1}{4}u^2T_cY^{(3,2)}(0,0))h_{\mathrm{}}^{\frac{8}{5}}+\mathrm{}\hfill \end{array}$$
(98)
## Appendix C Example of numerical data

Here we report, for sake of completeness, an example of the data we used to fit.
We refer to value of $`\beta =0.4373`$.


format: text

text:
# Evidence that an Intrinsic Redshift Component that is a Harmonic of z = 0.062 May be Present in Every Quasar Redshift

## 1 Introduction

It has been shown (Bell , 2002a, b) from the redshifts and positions of 12 of the compact, high-redshift objects near NGC 1068 that they appear to have been ejected outwards along the rotation axis of NGC 1068 in four similarly structured triplets that differ only in their size, rotation angle and rotation axis tip angle. In the previous studies (Bell , 2002a, b) no attempt was made to remove the redshift components due to the triplet ejection velocities. Assuming here that the model is correct, it is shown that it is possible to determine the rotation axis tip angle for each triplet at the time of its ejection. All ejection-related Doppler components can then be removed allowing the residual intrinsic redshift components to be obtained. This is the first instance where an attempt has been made to derive accurate values for the intrinsic components. Previously, they have always been contaminated with cosmological and other Doppler-related components. It is shown further that the intrinsic redshift values derived for the sources near NGC 1068 can be used to define a general relation, and that this relation, in turn, predicts the values found to be preferred in emission line redshift distributions covering redshifts from z = 0.06 to z = 2.

## 2 Calculating Intrinsic Redshift Components

The compact objects near NGC 1068 (Burbidge , 1999) are plotted in Fig. 1, where the heavy solid lines connect the single and paired source in each triplet, and the rotation axis of the central torus is indicated by the vertical dashed line. The four triplets are defined by the four ellipses labeled A to D. In the ejection model presented earlier (Bell , 2002a, b), each triplet is ejected outward along the rotation axis of the central torus. The axis is assumed to have an associated tip angle $`\gamma `$ defined as its angle of tip toward the observer at the top. The triplet center of mass (CofM) moves outward along this axis and its speed and direction must be known if the Doppler components resulting from this motion are to be calculated and removed. Because the tip angle was unknown, no attempt was made to remove this ejection component in the previous studies. In addition to this motion, the triplets further divide into a single and a pair that separate with equal velocities in opposing directions, each perpendicular to the direction of triplet ejection. The position angle defined by the line joining the single and pair-midpoint (the heavy lines in Fig. 1) gives the rotation angle ($`\theta `$) of the relevant triplet (where $`\theta `$ is measured counter-clockwise from North). It was shown previously (Bell , 2002a) how $`\theta `$ can be calculated if the rotation axis tip angle is known.

a) Estimating the tip angle $`\gamma `$.

Previously the tip angle for all four triplets was assumed to be $`\gamma =18\mathrm{°}`$. However, there appears to be strong evidence to indicate that it may be significantly larger than this for all triplets. Since the orbital plane is perpendicular to the rotation axis, if the triplet CofM is ejected along the rotation axis, then it is likely that the triplet plane will remain perpendicular to the direction of ejection. The plane defined by the pair and singlet in each triplet thus defines the orientation of the rotation axis (the tip angle) at the time the triplet was ejected.

In this scenario, if $`\gamma `$ = 0 for a given triplet, its three sources would all lie along an E-W (horizontal) line in Fig. 1. In triplets A and B, where $`\theta `$ is close to -90$`\mathrm{°}`$, any increase in $`\gamma `$ above 0 introduces a N-S separation on the sky in the relevant paired sources. For triplets A and B the tip angles are given by the angle $`\beta `$, defined as the angle between the horizontal and the line joining the two sources making up the pair. In triplet D, where $`\theta `$ is close to 0$`\mathrm{°}`$, an increase in $`\gamma `$ would introduce a N-S displacement in the relevant pair-single separation. Thus, for triplet D the tip angle is given approximately by the angle between the horizontal and the line joining the single and pair midpoint. In the case of triplet C, whose rotation angle lies closer to -45$`\mathrm{°}`$, the above approach is less accurate.

The $`\gamma `$ values obtained using the above technique for triplets A, B, C and D are listed in col 7 of Table 1. These values are all significantly larger than the value $`\gamma =18\mathrm{°}`$ assumed previously, and in closer agreement with the tip angle near $`60\mathrm{°}`$ currently measured for the optical axis of NGC 1068 (Schmitt and Kinney , 2000; Bock et al. , 2000). The tip angle for triplet A can be either $`55\mathrm{°}`$ or $`125\mathrm{°}`$. The correct value appears to be $`125\mathrm{°}`$. The ejection model requires this when the relative orientations of the red- and blue-shifted components making up the various pairs are taken into account. This conclusion also means that triplet A will have been ejected in the same jet direction as the other three triplets and not in the counter-jet direction as previously assumed. The CofM of triplet A will then move toward us instead of away, and this must be taken into account when the Doppler components are removed.

b) Calculating Doppler Components.

It is assumed here that the tip angles calculated above, which have large uncertainties ($`5\mathrm{°}10\mathrm{°}`$), can at least act as a rough indication of the true tip angle. It is further assumed that other physical relations, which must be satisfied by the ejection model, can be used to define the tip angle more accurately. For example, if the model is valid, the event ages derived should vary smoothly and, if the mass is constant, linearly with the rotation angle.

In order to examine how changes in the tip angle affect the age of each triplet, a spreadsheet was prepared that calculated new ages and intrinsic redshifts with the insertion of each new value of $`\gamma `$. Since the redshifts of sources 10 and 12 have not been measured it was also necessary to examine a range of redshifts for these sources. This too was accommodated in the spreadsheet calculations. Initial values were obtained by assuming that the pair-single ejection velocities of triplets C and D were similar to those measured for triplets A and B. Column 9 of Table 1 list the redshifts derived for sources 10 and 12. Table 2 lists various stages in the calculation of triplet ages and Table 3 summarizes the stages involved in calculating intrinsic redshifts.

The $`\gamma `$-values listed in column 2 of Table 2 and the redshifts listed for the C and D triplets in column 5 were arrived at using the values in col 7 of Table 1 as a guide and requiring the $`\theta `$ vs Age relation to be linear.

c) Rotation period

The resulting $`\theta `$ vs Age relation is shown in Fig. 2. The slope of this line gives a rotation period for the nucleus of NGC 1068 of $`1.8\times 10^6`$ yr, which is a factor of 3.3 shorter than that found by Alloin et al.,  (2001) for the central torus. This value is assumed to be more accurate than that given earlier (Bell , 2002a) since the earlier age calculations did not include accurate corrections for all ejection velocities. Since the age is calculated from the present position of the objects and their ejection velocity perpendicular to the l-o-s, this would have affected the previous age calculations.

## 3 Intrinsic Redshift Relation

The intrinsic redshifts calculated in Table 3 are plotted versus level (triplet) in Fig. 3 where the redshift spacing between levels can be seen to increase in multiples of the redshift increment $`\mathrm{\Delta }z`$ = 0.062. Before the Doppler components due to triplet ejection were removed here, the redshifts (zmean) were found to exhibit similar characteristics (Bell , 2002a), except in that case the redshift increments between levels were found to be multiples of $`\mathrm{\Delta }z`$ = 0.05. As long as the ejection velocities are small, and the rotation axis is tipped toward the observer, the slightly larger increment of $`\mathrm{\Delta }z`$ = 0.062 found here after removal of the Doppler components would be anticipated. It is therefore not surprising that a similar relation was visible in the raw data. This also reinforces the likelihood that the results obtained here have not been introduced by the above calculations. The most interesting new result is that when all Doppler components are accounted for, the intrinsic redshifts are all multiples of the redshift increment $`\mathrm{\Delta }`$ z = 0.062 first found by Burbidge over thirty years ago (Burbidge , 1968; Burbidge and Hewitt , 1990).

The intrinsic redshift values found for the four triplets near NGC 1068 are labeled A, B, C, and D in Fig. 3. Two additional intrinsic redshift values (0.62 and 0.31) can be predicted from Fig. 3 corresponding to levels 4 and 5. These appear to be the only other acceptable values since higher ones would be blueshifted.

These intrinsic redshift components can be defined by the relation:

z$`{}_{\mathrm{i}}{}^{}=0.062[10NM]`$ —————– (1)

where $`N`$ = 2, $`M_2=[n(n+1)/2]`$, and $`n`$ = 0, 1, 2, 3, 4, and 5.

For N = 2, this relation predicts intrinsic redshift values of 0.31, 0.62, 0.868, 1.054, 1.178, and 1.24, all separated by multiples of the redshift increment $`\mathrm{\Delta }z`$ = 0.062, and all occurring at harmonics of z = 0.062. Thus this increment is present in the intrinsic component for redshifts well above the z = 0.06 to 0.6 redshift range considered by Burbidge and Hewitt  (1990).

## 4 Preferred Emission-Line Redshift Values found in Previous QSO Redshift Analyses

Several investigators have carried out searches for periodicities in quasar emission-line redshift distributions. It is important to note that all of these were carried out using observed redshifts that can include large Doppler ejection and cosmological components in addition to an intrinsic component. More than thirty years ago Burbidge  (1968) pointed out that there were strong peaks near z = 0.061 and 1.95. Burbidge and Hewitt  (1990) showed that the z = 0.06 periodicity held up when a larger sample (z $``$ 0.6) was used. More extensive analyses of the raw redshift distributions have found other preferred redshift values. Karlsson  (1971, 1977); Burbidge  (1978); Burbidge and Napier  (2001) reported that the redshift distributions showed peaks near z = 0.061, 0.30, 0.60, 0.96, 1.41, 1.96, 2.36, 3.45 and 4.47. Karlsson  (1977) pointed out that these redshift peaks can be reasonably well represented by the relation $`\mathrm{\Delta }`$log(1+z) = 0.089.

The distribution found for QSO emission line redshifts below z = 1.24 (Karlsson , 1977) is reproduced in Fig. 4. This plot has been chosen because it represents one of the first investigations. Although several follow-up investigations have been reported (Arp et al. , 1990; Burbidge and Napier , 2001), the peaks below z = 2.0 remain relatively unchanged. These later investigations also represent cases where data have been selected after a relation has been previously found and this can sometimes lead to data choices that reinforce the original relation. No such bias can be present in the initial investigation.

In Fig. 4 the heavy vertical bars near the bottom of the figure indicate the locations of the peaks claimed to be periodic in log(1+z). Also included in the figure are vertical dotted lines indicating the locations of the intrinsic redshifts defined by eqn 1 ($`N`$ = 2) for the sources near NGC 1068 (hereafter referred to as zN=2 values). They range from the maximum value z = 1.24, down to z = 0.31 as $`n`$ varies from 0 to 5. All occur at pure harmonics of z = 0.062. It is clear from this figure that, even though the emission-line redshifts do not represent a purely intrinsic redshift distribution, there is excellent agreement between the peaks in the distribution and the intrinsic redshifts found from the QSOs near NGC 1068. In fact, for z $`1.24`$, the fit to the zN=2 redshifts is significantly better than the fit to the log(1+z) periodicity. This implies that the Doppler components present in the redshifts of the sources near each peak must be small. This, in turn, implies that they must be associated with reasonably nearby objects (z$`{}_{c}{}^{}<0.02`$), and have been expelled approximately perpendicular to the l-o-s. Those not expelled perpendicular to the l-o-s are assumed to be represented by those sources located below the curved baseline in Fig. 5.

In previous analyses (Karlsson , 1977; Arp et al. , 1990; Burbidge and Napier , 2001) it seems not to have been realized that the emission-line redshift peaks are located at harmonics of z = 0.062. Presumably this can be explained by the facts that the peaks are not periodic in z and that the redshifts used are observed values and not purely intrinsic ones, which may affect the peak positions slightly. However, it is worth noting also that the log(1+z) periodicity has nothing to support it, whereas the zN=2 redshifts have been predicted here using an analysis that is completely independent from the redshift distribution analysis. Furthermore, the log(1+z) relation has been known for 25 years and has not led to any other related discoveries. This and the fact that it is based on a distribution that can contain large Doppler components are arguments to question its appropriateness to correctly represent the intrinsic redshift relation.

## 5 Other Values of $`N`$

The objects clustered near z = 0.061 appear to be of a different class (many are Seyferts, N-gals, etc.) and are not included in the distribution in Fig. 4. On the other hand, the objects examined by Burbidge and Hewitt  (1990), which had redshifts up to z = 0.6, are distributed at single intervals of $`\mathrm{\Delta }`$z = 0.062. One might then speculate that these objects (with a maximum redshift of z = 0.62) represent the $`N`$ = 1 state of eqn 1. If so, for zN=1, from the results of Burbidge and Hewitt  (1990), it can be concluded that $`M_1=n`$ where $`n`$ can have the ten values 0,1,2….9. Eqn 1 can then be expressed in the form z$`{}_{\mathrm{i}}{}^{}=[0.620.062n`$\].

In Fig. 5 all 574 objects examined by Karlsson  (1977) have been plotted. The objects with measured redshifts above z = 1.24 may simply be zN=2 sources with large Doppler components that have yet to be removed, or, they may represent the $`N`$ = 3 state (see discussion below on gravitational redshifts for further comments on redshifts above z = 2).

We have thus seen from the results of Burbidge and Hewitt  (1990), that for the $`N=1`$ state, $`M_1`$ = n. From the NGC 1068 results for the $`N=2`$ state, $`M_2`$ = $`n(n+1)/2`$. Thus as $`N`$ changes from $`N`$ = 1 to $`N`$ = 2, $`M`$ goes from $`M_1=n`$ to $`M_2=n(n+1)/2`$. If we assume that the $`N`$ = 3 state behavior can be derived from the $`N`$ = 2 state in the same manner that the $`N`$ = 2 state came from the $`N`$ = 1 state, then for $`N`$ = 3, M$`{}_{3}{}^{}=[p(p+1)/2]`$, where p = $`[n(n+1)/2]`$. It can easily be shown that acceptable values of $`M`$ are then $`M_3`$ = 0,1,6, and 21 where $`n`$ = 0,1,2, and 3. This gives intrinsic redshift values for the $`N`$ = 3 state of 1.86, 1.798, 1.488, and 0.558. As can be seen in Fig. 5 these values show excellent agreement with all peaks below z = 2 that were not previously fitted by the zN=2 lines. This can be seen clearly in Fig. 6 where the curved baseline in Fig. 5 has been removed. The curved baseline, as postulated above, is assumed to represent those sources in which l-o-s Doppler components are sufficiently large to smear out the periodicity in the intrinsic components.

Eqn 1 can thus define all three $`N`$-states provided the appropriate values for $`M`$ are used in each case. For $`N`$ = 1, $`M_1`$ = 0,1,2,3,…9; For $`N`$ = 2, $`M_2`$ = 0,1,3,6,10,15; for $`N`$ = 3, $`M_3`$ = 0,1,6, and 21. These values are obtained directly from the relations $`M_1`$ = n; $`M_2=[n(n+1)/2]`$; $`M_3=[p(p+1)/2]`$ (where p = $`n(n+1)/2)`$. In all cases $`n`$ = 0, 1,2,3…,n. The redshifts predicted by these values of $`N`$ and $`M`$ all occur at harmonics of z = 0.062 and they show the same distribution in redshift as do the emission-line peaks.

## 6 Intrinsic Redshift vs $`M`$

In Fig. 7 the intrinsic redshift vs $`M`$ relation is compared for all three N-states. As can be seen, and can be predicted from eqn. 1, the intrinsic redshifts fall off vs $`M`$ with the slope -0.062 for all three states. This slope is thus given by the periodicity first found by Burbidge  (1968).

## 7 An Alternate Way of Expressing the Intrinsic Redshift Relation

Although eqn 1 may represent the simplest form of the intrinsic redshift relation, it is worth pointing out that it can also be expressed in the form zi = \[0.62$`N`$ \+ 0.062$`M`$\]. Here the constant $`10\times 0.062`$ has been expressed as 0.62. This number is identical to the theoretical maximum gravitational redshift (zgmax) obtained from the Schwartzschild interior solution for a perfect fluid sphere (Burbidge and Burbidge , 1967; Bondi , 1964; Buchdahl , 1959, 1966), and may suggest that the intrinsic redshift component is somehow produced gravitationally. This suggestion is also reinforced by the fact that there appears to be no correlation between the peaks reported above z = 2 and intrinsic redshifts that can be predicted for $`N`$-values above 3. If the intrinsic redshifts are of gravitational origin this would be consistent with the well-known z = 2 upper limit for gravitational redshifts (Burbidge and Burbidge , 1967). It is therefore tentatively assumed here that observed redshifts above z = 2 contain large Doppler components but their intrinsic components remain below zi = 1.86.

## 8 The Significance of the parameters $`N`$ and $`M`$

At the moment one can only speculate on the significance of the parameters $`N`$ and $`M`$. Since $`N`$ is constant for a given source, its value may be related to some parameter (e.g. mass) of the parent object or its nucleus. If so, this would seem to imply that the nuclear mass is quantized with at least three acceptable levels. Since the ejected objects all have similar masses initially and appear to expel mass (or energy) with time from the central QSO into a surrounding host, the discrete nature of the parameter $`M`$ may then suggest that this mass is ejected in discrete ”packets”.

As far as the formation of galaxies is concerned, these results do not change significantly the picture discussed earlier (Bell , 2002a, b). They still suggest that galaxies are created continuously throughout the life of the Universe. However, if the intrinsic component is gravitationally produced, it would seem less likely that the model of Narlikar and Das  (1980) plays any role. In Fig. 8 the zN=2 intrinsic redshifts are plotted vs Age on a logarithmic scale. It is clear that these intrinsic components play a very short-lived part in the life of a galaxy, lasting only for the first 3 million years.

## 9 Conclusions

By removing all ejection-related Doppler components from the redshifts of the compact QSOs near NGC 1068 it has been shown that the remaining intrinsic components can be defined by the relation zi = $`0.062[10NM]`$ where $`N`$ = 2, $`M=[n(n+1)/2]`$ and $`n`$ = 0,1,2,3,4,5. This is the first instance in which it has been possible to separate out the intrinsic component. It is shown that this same relation will fit the results obtained by Burbidge and Hewitt  (1990) for redshifts between z = 0.06 and z = 0.6, if $`N`$ = 1 and $`M=n`$, where $`n`$ = 0,1,2,… 9. The $`M`$-values for the $`N`$ = 3 state follow directly from these two states and it is shown that together these three states predict intrinsic redshifts that agree closely with those values found to be preferred in emission-line redshift analyses from z = 0.06 to z = 2. These results indicate that the $`\mathrm{\Delta }`$z = 0.062 redshift interval may thus represent a fundamental redshift spacing that arises from the very nature of the intrinsic redshifts.

Because this relation includes a constant (zgmax = 0.62) which can be related to gravitational redshifts, the intrinsic components may have a gravitational origin. If so, this increment may be evidence for some form of gravitational quantization.

It is concluded here that the log(1+z) relation, previously fitted to the peaks in the emission-line redshift distribution, may not be the appropriate one to describe the intrinsic redshifts.

By fine-tuning the tip angles it has also been possible to predict accurate redshift values for sources 10 and 12. It is important to note that values of the Doppler components calculated here could change if the redshifts of sources 10 and 12 can be measured and turn out to be significantly different from those estimated above. Small differences might still be accommodated by adjusting the relevant tip angle, provided the age versus rotation angle relation in Fig. 2 is not destroyed. However, the excellent agreement between the intrinsic redshifts given by eqn 1, and the redshift peaks found in quasar emission-line analyses, strongly suggests that the above relation is a significant one.