text:
# Classical stability of U(1)A domain walls in dense matter QCD

## I Introduction

In general, there do not exist domain walls or other topological defects within the Standard Model. This is due to the trivial topology of the vacuum manifold. In contrast, these objects are quite common in condensed matter physics and cosmology. However, it has been realized only recently  that topological defects such as domain walls and strings may exist within the Standard Model at a large chemical potential.

It is well known that topological defects result from symmetries being broken. In the past few years there has been renewed interest in high density QCD. Similar to the BCS pairing in conventional superconductivity, the ground state of QCD at high density is unstable due to the formation of a diquark condensate  (see  for a review). In this new ground state, various symmetries which are present at $`\mu =0`$ are broken by the presence of this non-zero diquark condensate. This leads to the formation of the various topological defects discussed in .

In , it was shown that at high densities ($`\mu 1\mathrm{GeV}`$) there exist domain wall solutions which interpolate between the same vacuum state, in which the U(1)A phase of the diquark condensate varies between $`0`$ and $`2\pi `$. This type of domain wall which interpolates between the same vacuum state has been studied before in the context of axion models . It is interesting to note that similar domain wall configurations are present for the zero density case in the large $`N_c`$ limit and could be present for the physically relevant case of $`N_c=3`$ , in which case they can be studied at the BNL Relativistic Heavy Ion Collider (RHIC) . Given this, one might ask what happens between these two regions of $`\mu =0`$ and $`\mu 1\mathrm{GeV}`$. The main goal of the present paper is the analysis of the classical stability of the U(1)A domain walls for $`\mu 1\mathrm{GeV}`$. We demonstrate that the U(1)A domain walls are classically stable down to densities $`\mu 800\mathrm{MeV}`$. The main idea is to interpolate between $`\mu 1\mathrm{GeV}`$ (where perturbation theory is justified) and $`\mu \mu _c`$ (where instanton calculations lead to a reasonable description). Ideally, we would like to be able to make definitive statements on the stability of such configurations all the way down from large $`\mu `$ to $`\mu _c500\mathrm{MeV}`$, where $`\mu _c`$ is the critical chemical potential above which the color superconducting phase occurs. Unfortunately, we are unable to make a definitive statement for $`\mu _c\mu 800\mathrm{MeV}`$ due to a lack of theoretical control in this region. One can only speculate on the behavior in this region.

This work is organized as follows. In Sec. 2, we will construct the effective potential. In Sec. 3, we will describe the nontrivial domain wall solutions. In Sec. 4, we will examine the classical stability of such configurations under small perturbations. In Sec. 5, we will end with concluding remarks and future considerations.

## II The Effective Potential

It is well known that the ground state of $`N_f=2,3`$, $`N_c=3`$ QCD exhibits the Cooper pairing phenomenon as in conventional superconductivity (known as the color 2SC (2 flavor color superconducting) and CFL (color-flavor locked) phases of QCD ). In what follows, we will consider the $`N_c=N_f=3`$ CFL phase. In the CFL phase, the condensates take the form

$`q_{L\alpha }^{ia}q_{L\beta }^{jb}^{}`$ $``$ $`ϵ_{\alpha \beta \gamma }ϵ^{ij}ϵ^{abc}X_c^\gamma ,`$ (1)
$`q_{R\alpha }^{ia}q_{R\beta }^{jb}^{}`$ $``$ $`ϵ_{\alpha \beta \gamma }ϵ^{ij}ϵ^{abc}Y_c^\gamma ,`$ (2)
where $`L`$ and $`R`$ represent left and right handed quarks; $`\alpha `$, $`\beta `$, and $`\gamma `$ are the flavor indices; i and j are spinor indices; a, b, and c are color indices; and $`X_c^\gamma `$ and $`Y_c^\gamma `$ are the condensates which are complex color-flavor matrices. As the magnitudes of these condensates depend on the color index $`c`$, one can easily see that these objects are not gauge invariant by themselves. In order to construct a gauge invariant field, the following matrix which describes the octet of mesons and axial singlet was considered in :

$$\mathrm{\Sigma }_\gamma ^\beta =XY^{}=\underset{c}{}X_c^\beta Y_\gamma ^c.$$
(3)
If a U(1)A rotation ($`q\mathrm{exp}^{i\gamma _5\alpha /2}q`$) of this gauge-invariant field $`\mathrm{\Sigma }`$ is performed, we see that the fields (1) transform as

$`Xe^{i\alpha }X,`$ (4)
$`Ye^{+i\alpha }Y,`$ (5)
and therefore

$$\mathrm{\Sigma }e^{2i\alpha }\mathrm{\Sigma }.$$
(6)
Goldstone’s theorem states that there must be a single Goldstone mode $`\eta ^{}`$ associated with the breaking of this symmetry. Given this, we can parametrize the field as follows:

$$\mathrm{\Sigma }=\mathrm{\Sigma }_oe^\rho e^{i\varphi },$$
(7)
where the phase $`\varphi =\eta ^{}/f`$ is defined as a dimensionless field which describes the $`\eta ^{}`$ boson, $`f`$ is the corresponding decay constant, $`\mathrm{\Sigma }_o`$ is the vacuum expectation value of the composite field (3), and $`\rho `$ is another dimensionless field which describes the fluctuations of the magnitude of the condensate***The field $`\rho `$ should not be confused with the familiar $`\rho `$ meson in QCD. (analogous to the $`\sigma `$ field related to the fluctuations of the $`\overline{q}q`$ chiral condensate). We choose to parametrize the field as $`e^\rho `$ for convenience later on.

In order to construct an effective potential describing the dynamics of the phase of the condensate $`\varphi `$ as well as the magnitude $`|\mathrm{\Sigma }|`$, there are two types of terms which must be included. The first term which explicitly breaks the U(1)A symmetry was calculated in  by substituting the form of condensates given above into the instanton induced four-fermion Lagrangian :

$$V_{1\mathrm{inst}}(\rho ,\varphi )=a\mu ^2\mathrm{\Delta }_o^2e^\rho \mathrm{cos}\varphi ,$$
(8)
where $`\mathrm{\Delta }_o`$ is the value of the gap in the quark spectrum. The perturbative form of the expectation value of the condensate has been used in arriving at this result :

$$|\mathrm{\Sigma }|\mathrm{\Sigma }_o=\frac{9}{8\pi ^2}\frac{\mu ^4\mathrm{\Delta }_o^2}{g^2}.$$
(9)
For $`N_f=3`$, the dimensionless coefficient $`a`$ was found to be

$$a=7\times 10^3\left(\frac{m_s}{\mu }\right)\left(\mathrm{ln}\frac{\mu }{\mathrm{\Lambda }_{\mathrm{QCD}}}\right)^7\left(\frac{\mathrm{\Lambda }_{\mathrm{QCD}}}{\mu }\right)^9,$$
(10)
where $`m_s`$ is the mass of the strange quark. The mass of the corresponding $`\eta ^{}`$ boson can be easily calculated by expanding the potential, $`m_\eta ^{}=2\pi \sqrt{a}\mathrm{\Delta }_o`$. According to , this potential is only under theoretical control when the mass of the $`\eta ^{}`$ boson is much less than the typical scale for higher excitations $`2\mathrm{\Delta }`$, which corresponds to $`a1/\pi ^2`$. For physical values of the strange quark, this corresponds to a chemical potential of about $`\mu 700\mathrm{MeV}`$ (for $`a1/\pi ^2`$).

To be able to draw any conclusions about the stability of the domain wall solution for intermediate densities $`(\mu _c\mu 1\mathrm{GeV})`$, one must include degrees of freedom which are related to the fluctuations of the absolute value of $`|\mathrm{\Sigma }|`$. We do not know the effective potential in the region of interest; however, for qualitative discussions we shall use a potential derived for asymptotically large $`\mu `$. The second type of term which must be included in an effective potential description is one which uniquely fixes the magnitude of the condensate. An effective potential for the magnitude of the condensate $`|\mathrm{\Sigma }|`$ was derived in  in the perturbative region where the analytical form of the gap is known. This effective potential fixes uniquely the value of the vacuum condensate. The potential is of the Coleman-Weinberg type  and is given as follows:

$`V_{\mathrm{pert}}(|\mathrm{\Sigma }|)`$ $`=`$ $`\left({\displaystyle \frac{2\nu \pi }{\mu }}\right)^2|\mathrm{\Sigma }|\left(1\mathrm{ln}{\displaystyle \frac{|\mathrm{\Sigma }|}{\mathrm{\Sigma }_o}}\right),`$ (11)
$`=`$ $`{\displaystyle \frac{\mu ^2\mathrm{\Delta }_o^2}{\pi ^2}}e^\rho (1\rho ),`$ (12)
$`\nu `$ $`=`$ $`\sqrt{{\displaystyle \frac{8\alpha _s}{9\pi }}},`$ (13)
where $`\alpha _s`$ is the standard strong coupling constant and the perturbative result of the condensate (9) has been used. In all calculations that follow we will assume $`\mathrm{\Delta }_o=100\mathrm{MeV}`$ as the numerical value for the gap. It should be noted that the potential (11) is justified only in the region

$$\nu \mathrm{log}(|\mathrm{\Sigma }|/\mathrm{\Sigma }_o)1.$$
(14)
Since we are considering $`|\mathrm{\Sigma }|\mathrm{\Sigma }_o`$, the use of this effective potential is justified.

We are interested in the region $`\mu _c\mu 1\mathrm{GeV}`$, where Eq. (11) is not literally correct. Even though the results stated above are not necessarily under theoretical control for $`\mu \mu _c`$, one can speculate on how the coefficients behave at intermediate densities when $`\mu _c\mu 1\mathrm{GeV}`$. Due to the fact that all the same symmetries are present as the chemical potential is lowered until reaching $`\mu _c`$, we expect the qualitative form of the effective potential (8) and (11) to remain the same. As the chemical potential is lowered, eventually perturbative calculations which are valid at asymptotically large $`\mu `$ are no longer the correct description and instanton calculations become relevant. One would expect that these calculations must match up at some point. However, the coefficients in front of the potential could possibly be very different in the density region of interest. We will refer to the coefficient in front of the one-instanton potential as $`\beta _1`$ and the coefficient in front of the perturbative potential as $`\beta _2`$. The value of $`\beta _1`$ is essentially fixed by the form of the condensate and by the constituent quark mass. Below the critical chemical potential $`\mu _c`$ at which the chiral phase transition occurs, this coefficient is independent of $`\mu `$ . Therefore, we would expect that $`\beta _1`$ would be some smooth function of $`\mu `$ which reaches its maximum value at $`\mu =\mu _c`$. At $`\mu _c\mu 1\mathrm{GeV}`$, the coefficient $`\beta _2`$ is modified by the formation of instanton-antiinstanton ($`I\overline{I}`$) molecules . We should also note that the coefficient $`\beta _1`$ can also be estimated from the instanton liquid model  using the average size of instantons as well as requiring a constituent quark mass of about $`350400\mathrm{MeV}`$.

Combining both terms we have an effective potential which is given by the following:

$$V(\rho ,\varphi )=\beta _2e^\rho (1\rho )\beta _1e^\rho \mathrm{cos}\varphi .$$
(15)
The values of the coefficients $`\beta _1`$ and $`\beta _2`$ are known for asymptotically large $`\mu `$:

$`\beta _1`$ $`=`$ $`a\mu ^2\mathrm{\Delta }_o^2,`$ (16)
$`\beta _2`$ $`=`$ $`{\displaystyle \frac{\mu ^2\mathrm{\Delta }_o^2}{\pi ^2}}.`$ (17)
In Fig. 1, we show the cross section of the potential at $`\varphi =0,2\pi `$ and $`\varphi =\pi `$ for $`\mu =800\mathrm{MeV}`$. Notice that the existence of an absolute minima at $`\varphi =0,2\pi `$ and a saddle point at $`\varphi =\pi `$ allows for nontrivial configurations which wind around the barrier at $`|\mathrm{\Sigma }|=0`$. In the limit $`\mu \mathrm{}`$, the parameter $`a0`$ and the potential has degenerate minima at $`|\mathrm{\Sigma }|=\mathrm{\Sigma }_o`$. The kinetic term is given by

$$\frac{|_o\mathrm{\Sigma }|^2u^2|_i\mathrm{\Sigma }|^2}{|\mathrm{\Sigma }|^2}=(_o\rho )^2u^2(_i\rho )^2+(_o\varphi )^2u^2(_i\varphi )^2,$$
(18)
where $`u`$ is the velocity which is different from $`1`$. The perturbative values for the decay constant $`f`$ and velocity $`u`$ were calculated in . In order to fix the correct dimensionality, we must multiply the kinetic term by the appropriate powers of the decay constant. The full effective Lagrangian up to two derivatives in the fields is then given by

$$=f^2[(_o\rho )^2u^2(_i\rho )^2]+f^2[(_o\varphi )^2u^2(_i\varphi )^2]V(\rho ,\varphi ).$$
(19)
In the above we have assumed that $`f_\eta ^{}f_\rho =f`$. The exact numerical value for $`f_\rho `$ is not known. However, in the large $`\mu `$ limit $`f_\eta ^{}f_\rho \mu `$ in order to have an appropriate scale for $`m_\rho \mathrm{\Delta }_o`$ (once again, the $`\rho `$ field should not be confused with the well known $`\rho `$ meson in QCD at $`\mu =0`$).

## III Domain wall solutions

As was done in , if we replace the field $`|\mathrm{\Sigma }|=\mathrm{\Sigma }_oe^\rho `$ by its vacuum expectation value (which is justified for $`\mu 1\mathrm{GeV}`$), the resulting potential is of the sine-Gordon type. The Lagrangian is given by the following:

$$=f^2[(_0\varphi )^2u^2(_i\varphi )^2]V_{1\mathrm{inst}}(\varphi ),$$
(20)
where the constant term has been dropped. The static domain wall solution to the corresponding equation of motion is well known. Considering a domain wall in the $`z`$ direction, the solution is given by

$$\varphi _o(z)=4\mathrm{arctan}\left(\mathrm{exp}(mz/u)\right),$$
(21)
where $`m`$ is the mass of the $`\eta ^{}`$. This solution interpolates between the same vacuum state; at $`z=\pm \mathrm{}`$ we have $`\varphi =0,2\pi `$, respectively. It is well known that this solution is absolutely stable under small perturbations $`\varphi =\varphi _o+\delta \varphi `$. In other words, the Schrödinger-type equation obtained by varying the field and linearizing the equation of motion,

$$_z^2\delta \varphi _n+\frac{m^2}{u^2}\left(12\mathrm{s}ech^2\left(\frac{m}{u}z\right)\right)\delta \varphi _n=\omega _n^2\delta \varphi _n,$$
(22)
has the lowest eigenvalue $`\omega _o=0`$ corresponding to $`\delta \varphi _o=d\varphi _o(z)/dz\mathrm{s}ech(mz/u)`$. This is just the zero mode which is a result of translational invariance $`zz+z_o`$. Since the lowest eigenvalue is non-negative, the domain wall solution is stable under small perturbations. It turns out that this is the only bound state which is a solution of Eq. (22).

In the case in which the replacement $`|\mathrm{\Sigma }|\mathrm{\Sigma }_o`$ is not done, the solution must be modified. If we want to study a stable solution for $`\mu `$ which is not asymptotically large, we must include fluctuations in $`\varphi `$ as well as in the $`\rho `$ direction (i.e., the absolute value of $`\mathrm{\Sigma }`$). The two equations of motion for static solutions are given by

$`2f^2u^2_i^2\rho `$ $`=`$ $`\beta _2\rho e^\rho \beta _1e^\rho \mathrm{cos}\varphi ,`$ (23)
$`2f^2u^2_i^2\varphi `$ $`=`$ $`\beta _1e^\rho \mathrm{sin}\varphi .`$ (24)
Although we do not know the exact solution for this set of coupled nonlinear differential equations, if $`\beta _1/\beta _2<1`$ we can approximate the solutions. In this case, the approximate solutions can be parametrized by

$`e^{\rho _o}`$ $``$ $`1+\alpha \mathrm{cos}\varphi _o,`$ (25)
$`\varphi _o`$ $``$ $`4\mathrm{arctan}\left(\mathrm{exp}(mz/u)\right),`$ (26)
where $`\alpha \beta _1/\beta _2`$. Our stability analysis will be based upon these approximate solutions of the equations of motion. Since the above solutions do not correspond to the exact solutions to the equations of motion (minimum energy path which winds around that barrier at $`|\mathrm{\Sigma }|=0`$), there will be nonzero linear terms when the energy of the system is perturbed about the domain wall solutions $`\varphi _o`$ and $`\rho _o`$. These will be estimated in the following section where the stability analysis is performed.

## IV Stability analysis of domain walls

Although these domain walls may exist as classically stable objects at large densities, it is not immediately obvious if this is the case at intermediate densities. In order to examine the classical stability of the domain wall configurations in the region $`\mu 1\mathrm{GeV}`$, we must look at how the system reacts to small perturbations of the fields. The energy density is given by the following expression:

$$E(\varphi ,\rho )=_{\mathrm{}}^+\mathrm{}𝑑z[f^2u^2(\rho )^2+f^2u^2(\varphi )^2+V(\varphi ,\rho )].$$
(27)
In order to express all integrals as dimensionless quantities, we will perform the following change of variables:

$$z^{}=\frac{z}{\lambda },\lambda =\frac{u}{m_\eta ^{}}.$$
(28)
The energy density is now given by

$$E(\varphi ,\rho )=\frac{f^2u^2}{\lambda }\left[𝑑z^{}\left((\rho )^2+(\varphi )^2+\frac{\lambda ^2}{f^2u^2}V(\varphi ,\rho )\right)\right],$$
(29)
where the derivative is now taken with respect to the dimensionless coordinate $`z^{}`$. We can see that this has the correct dimensions of $`\mathrm{MeV}^3`$ times the dimensionless integral in square brackets. Following the standard method for analyzing the stability of a classical solution which was briefly described in the previous section, we will expand the fields about their vacuum expectation values:

$`\rho \rho _o+\delta \rho ,`$ (30)
$`\varphi \varphi _o+\delta \varphi .`$ (31)
Next, we substitute Eq. (30) into Eq. (27) and perform an expansion, neglecting any terms greater that quadratic order in $`\delta \rho `$ and $`\delta \varphi `$:

$`E`$ $``$ $`E^{(0)}+E^{(1)}+E^{(2)}`$ (32)
$``$ $`E(\varphi _o,\rho _o)+\gamma {\displaystyle 𝑑z^{}\delta \rho \left(2^2\rho _o+\frac{\lambda ^2}{f^2u^2}\frac{\delta V}{\delta \rho }|_{\varphi _o,\rho _o}\right)}`$ (33)
$`+`$ $`\gamma {\displaystyle 𝑑z^{}\delta \varphi \left(2^2\varphi _o+\frac{\lambda ^2}{f^2u^2}\frac{\delta V}{\delta \varphi }|_{\varphi _o,\rho _o}\right)}+\gamma {\displaystyle 𝑑z^{}\delta \rho \left(^2+\frac{\lambda ^2}{2f^2u^2}\frac{\delta ^2V}{\delta \rho ^2}|_{\varphi _o,\rho _o}\right)\delta \rho }`$ (34)
$`+`$ $`\gamma {\displaystyle 𝑑z^{}\delta \varphi \left(\frac{\lambda ^2}{f^2u^2}\frac{\delta ^2V}{\delta \varphi \delta \rho }|_{\varphi _o,\rho _o}\right)\delta \rho }+\gamma {\displaystyle 𝑑z^{}\delta \varphi \left(^2+\frac{\lambda ^2}{2f^2u^2}\frac{\delta ^2V}{\delta \varphi ^2}|_{\varphi _o,\rho _o}\right)\delta \varphi },`$ (35)
where $`\gamma =f^2u^2/\lambda `$. The first term $`E^{(0)}`$ in the above expansion is the energy density or wall tension of the domain wall. In the case in which the domain wall solutions given by Eq. (25) were the exact solutions to the classical equations of motion, the linear terms $`E^{(1)}`$ (proportional to $`\delta \rho `$ and $`\delta \varphi `$) would be zero everywhere. Due to the fact that our solutions are not exact, these must be considered.

First, let us estimate the term which is linear in $`\delta \varphi `$:

$$\delta \varphi ^{(1)}=\frac{f^2u^2}{\lambda }𝑑z^{}\delta \varphi \left(2^2\varphi _o+\frac{\lambda ^2}{f^2u^2}\beta _1e^{\rho _o}\mathrm{sin}\varphi _o\right).$$
(36)
We know that $`e^{\rho _o}1+\rho _o`$ and using the fact that $`\varphi _o`$ is a solution to the equation of motion given by Eq. (24) with $`\rho =0`$, we have:

$`\delta \varphi ^{(1)}`$ $``$ $`{\displaystyle \frac{f^2u^2}{\lambda }}{\displaystyle 𝑑z^{}\frac{\lambda ^2}{f^2u^2}\beta _1\rho _o\mathrm{sin}\varphi _o\delta \varphi }`$ (37)
$``$ $`{\displaystyle \frac{f^2u^2}{\lambda }}{\displaystyle 𝑑z^{}a\pi ^2\mathrm{sin}2\varphi _o\delta \varphi }.`$ (38)
This linear term goes like $`a\pi ^2`$ ($`0.3`$ for $`\mu =800\mathrm{MeV}`$) and the integrand is small compared to the wall tension $`E^{(0)}`$ and can be neglected. The term which is linear in $`\delta \rho `$ is:

$$\delta \rho ^{(1)}=\frac{f^2u^2}{\lambda }𝑑z^{}\delta \rho \left(2^2\rho _o+\frac{\lambda ^2}{f^2u^2}\beta _2\rho _oe^{\rho _o}\frac{\lambda ^2}{f^2u^2}\beta _1e^{\rho _o}\mathrm{cos}\varphi _o\right).$$
(39)
Using again that $`\rho _o(\beta _1/\beta _2)\mathrm{cos}\varphi _o`$ we see this simplifies to

$`\delta \rho ^{(1)}`$ $``$ $`{\displaystyle \frac{f^2u^2}{\lambda }}{\displaystyle 𝑑z^{}\left(2a\pi ^2\delta \rho ^2(\mathrm{cos}\varphi _o)\right)}.`$ (40)
Since this term is also proportional to $`a\pi ^2`$, $`\delta \rho \rho _o`$, and the integral of $`^2(\mathrm{cos}\varphi _o)`$ is small, we can also neglect this term. The magnitude of the linear term shows how far away we are from the exact solution. This information will be used in what follows for the stability analysis.

Now we consider the most important quadratic term. A similar case involving two coupled scalar fields was looked at in  and we will follow the standard procedure presented there closely. If the field configuration is classically stable, the second variation of the energy should be a positive differential operator. This means we must solve the following Schrödinger-type eigenvalue problem:

$$H\left(\begin{array}{c}\delta \rho \\ \delta \varphi \end{array}\right)=\omega ^2\left(\begin{array}{c}\delta \rho \\ \delta \varphi \end{array}\right),$$
(41)
where $`H`$ is the operator,

$$H=_z^{}^2\mathrm{𝟏}+\frac{\lambda ^2}{2f^2u^2}U,$$
(42)
and 1 is the $`2\times 2`$ identity matrix. The potential $`U`$ is a $`2\times 2`$ matrix with elements:

$`U_{11}`$ $`=`$ $`\beta _2(1+\rho _o)e^{\rho _o}\beta _1e^{\rho _o}\mathrm{cos}\varphi _o,`$ (43)
$`U_{12}`$ $`=`$ $`U_{21}=\beta _1e^{\rho _o}\mathrm{sin}\varphi _o,`$ (44)
$`U_{22}`$ $`=`$ $`\beta _1e^{\rho _o}\mathrm{cos}\varphi _o.`$ (45)
If the domain wall solution is a stable one, then the operator $`H`$ is positive-semidefinite. The eigenvalue equations can be decoupled by diagonalizing the matrix $`U`$. We should note that only the potential term has to be diagonalized when looking for negative energy modes, as was done in . The result is

$`U_\pm ={\displaystyle \frac{1}{2}}\left(a\pm \sqrt{a^24(abb^2c^2)}\right),`$ (46)
where $`a`$, $`b`$, and $`c`$ are defined as

$`a`$ $`=`$ $`\beta _2(1+\rho _o)e^{\rho _o},`$ (47)
$`b`$ $`=`$ $`\beta _1e^{\rho _o}\mathrm{cos}\varphi _o,`$ (48)
$`c`$ $`=`$ $`\beta _1e^{\rho _o}\mathrm{sin}\varphi _o.`$ (49)
The operator $`H`$ now takes the following form:

$$H=\left(\begin{array}{cc}_z^{}^2+\frac{\lambda ^2}{2f^2u^2}U_+& 0\\ 0& _z^{}^2+\frac{\lambda ^2}{2f^2u^2}U_{}\end{array}\right).$$
(50)
Since we can immediately see that $`U_+0`$ for all $`z^{}`$ due to the fact that $`a>0`$, there does not exist any negative eigenvalues corresponding to the first equation in this transformed basis. Looking at the second eigenvalue equation, we have

$$\left(_z^{}^2+\frac{\lambda ^2}{2f^2u^2}U_{}\right)\psi _n=\omega _n^2\psi _n.$$
(51)
It is a well known theorem of quantum mechanics that there must exist at least one bound-state solution to Eq. (51). Due to the fact that our domain wall solution should be invariant under translations in space $`z^{}z^{}+z_o^{}`$, there should be a corresponding zero mode in the spectrum of Eq. (51). In the high density limit, we recover the familiar sine-Gordon equation and we know that there is only one bound state in the spectrum of Eq. (22). If the exact solution to the equations of motion (23) and (24) were known, one would expect to see a corresponding mode with a vanishing eigenvalue in the spectrum of Eq. (51). As the density is lowered ($`\beta _1\beta _2`$) and the saddle point at $`\varphi =\pi `$ is still present, one would expect the appearance of a mode with a negative eigenvalue corresponding to instability of the domain wall. Due to the fact that there still must be a zero mode in the spectrum, the zero mode would become the first excited state of Eq. (51) and the lowest mode would have some negative eigenvalue $`\omega _0^2<0`$ corresponding to the instability of the domain wall. The problem of stability analysis now reduces to determining the eigenvalues corresponding to the bound states of Eq. (51). The appearance of an additional bound state in the spectrum as the chemical potential is lowered will be the first sign that the system is approaching the point of instability.

Since the solution corresponding to Eq. (25) is not the exact solution but does represent a path which winds around the barrier at $`|\mathrm{\Sigma }|=0`$, it is quite possible that the zero mode could show up in the spectrum with a small nonzero eigenvalue. It would show up as a true zero mode only when the exact solution to the equation of motion is substituted into Eqs. (47).

Although the potential $`U_{}`$ is nontrivial, we will use a variational approach in order to determine the upper bounds on $`\omega _0^2`$ and $`\omega _1^2`$. In choosing a trial wave function, we make the observation that the potential $`U_{}`$ is quite similar to the same potential which arises when analyzing the stability of the sine-Gordon soliton, Eq. (22). We will pick our normalized trial wave function accordingly:

$$\psi _0=\sqrt{\frac{\sigma }{2}}\mathrm{s}ech(\sigma z^{}),$$
(52)
with $`\sigma `$ being the variational parameter. Note that this trial wave function satisfies the required boundary conditions $`\psi _o(z^{}=\pm \mathrm{})0`$. For the first excited state, we must pick an odd function of $`z`$. We will choose the properly normalized function

$$\psi _1=\sqrt{\frac{3\sigma }{2}}\mathrm{tanh}(\sigma z^{})\mathrm{s}ech(\sigma z^{}).$$
(53)
Applying the variational principle, we must calculate

$$E_n^{(2)}(\sigma )=\frac{f^2u^2}{\lambda }\psi _n|\left(_z^{}^2+\frac{\lambda ^2}{2f^2u^2}U_{}\right)|\psi _n,$$
(54)
and minimize this quantity with respect to $`\sigma `$ to obtain an upper bound $`E_n^{(2)}(\sigma )`$ on the energy of the $`n\mathrm{t}h`$ state. The integral given by Eq. (54) must be performed numerically.

We now have the framework in place in order to test the stability of our domain wall configurations. We will assume the perturbative value for the constants $`f`$ and $`u`$ as calculated in , $`f^2=\mu ^2/(8\pi ^2)`$ and $`u^2=\frac{1}{3}`$. Setting $`\mu =800\mathrm{MeV}`$ and $`m_s=100\mathrm{MeV}`$, we see that the ratio of the coefficients is $`\beta _1/\beta _20.3`$. In this case, the linear terms (37) and (40) are small and the solutions given by Eq. (25) are valid approximations to the exact solutions. In Fig. 2, we show the effective Schrödinger potential $`U_{\mathrm{eff}}=\lambda ^2U_{}/(2f^2u^2)`$ for $`\mu =800\mathrm{MeV}`$. For the above choice of parameters, the wall tension given by Eq. (29) is

$$E^{(0)}(\varphi _o,\rho _o)17.48\frac{f^2u^2}{\lambda }.$$
(55)
For the trial wave functions given in Eqs. (52) and (53), the following results for the two bound states of Eq. (51) were obtained:

$`E_0^{(2)}(\sigma _{\mathrm{min}})0.016{\displaystyle \frac{f^2u^2}{\lambda }},`$ (56)
$`E_1^{(2)}(\sigma _{\mathrm{min}})+1.163{\displaystyle \frac{f^2u^2}{\lambda }}.`$ (57)
From this, we can see that $`E_1^{(2)}E_0^{(2)}`$ and both of these quantities are much less the wall tension given by Eq. (55). Even though the ground state energy seems to be negative, due to the fact that $`E_1^{(2)}E_0^{(2)}`$ we can associate this mode with the zero mode. The small nonzero eigenvalue is actually an artifact of our approximations. The appearance of a negative mode is merely a consequence of the approximate solutions (25) as discussed above. This identification can be verified by increasing the chemical potential. When the above calculations are repeated as the chemical potential is increased, we see that the energy $`E_0^{(2)}0`$. This result is expected due to the fact that the $`|\mathrm{\Sigma }|`$ field can be integrated out as $`\mu `$ increases and the sine-Gordon type theory is recovered. As $`\mu `$ is increased, we also see that the eigenvalue of the first excited state $`\omega _1`$ increases towards the maximum value of $`U_{\mathrm{eff}}`$. Eventually, as $`\mu `$ is increased further, there is no longer a first excited bound state in the spectrum. The calculation of $`E_0^{(2)}`$ and $`E_1^{(2)}`$ was done with the variational functions chosen to be different from Eqs. (52) and (53). The results obtained were the same order of magnitude as stated above (56). This supports our interpretation of the $`E_0^{(2)}`$ state as the would-be translational mode if exact solutions are known. The magnitudes of the linear terms (37) and (40) are approximately the same order of magnitude of $`E_0^{(2)}`$, which supports our interpretation of $`E_0^{(2)}`$ as the zero mode.

## V Conclusion

The main goal of this paper was an analysis of the classical stability of U(1)A domain walls . Naively one would expect that decreasing $`\mu `$ from $`\mu 1\mathrm{GeV}`$ (when the calculations are under control ), we inevitably face the situation in which the domain walls become unstable objects due to the fast growth of the coefficient $`\beta _1`$ (15). This naive expectation may not necessarily be correct due to the even faster growth of the coefficient $`\beta _2`$ (15), which receives contributions from the formation of $`I\overline{I}`$ molecules as well as perturbative contributions.

What we have actually demonstrated is that the domain wall solution remains classically stable down to $`\mu 800\mathrm{MeV}`$. In order to assess the stability of the domain walls for $`\mu _c\mu 800\mathrm{MeV}`$, one must explicitly include the $`I\overline{I}`$ contribution in the effective potential (which is expected to be the dominant contribution at $`\mu \mu _c`$ ). This would hopefully lead to $`\beta _1/\beta _2<1`$, ensuring the classical stability of the U(1)A domain walls. Unfortunately, due to the lack of information in this region we cannot generalize our results to below $`800\mathrm{MeV}`$ to $`\mu _c`$. We can argue that the U(1)A domain walls remain classically stable down to $`\mu \mu _c`$ due to the faster growth of the coefficient $`\beta _2`$ compared to $`\beta _1`$ as $`\mu `$ is decreased. The ratio $`\beta _1/\beta _2`$ may be very sensitive to changes in the instanton size distribution or the various form factors. This is a difficult problem and unfortunately I do not know how to estimate such contributions for a small chemical potential.

We should remark here that the stability of the U(1)A domain wall implies a classical stability of the U(1)A strings , which become the edge of the domain walls. It remains to be seen whether these (or other) topological defects will have any impact on the physics of neutron stars and other compact stellar objects with high core density $`\mu `$.

###### Acknowledgements.

I would like to thank A. Zhitnitsky for suggesting this problem to me and for many helpful discussions and suggestions. I would also like to thank T. Schäfer for interesting discussions during the National Nuclear Physics Summer School in Bar Harbor, ME. This work was supported in part by the Natural Sciences and Engineering Research Council of Canada.


format: text

text:
# Structure of the exotic spin-flop states in BaCu2Si2O7.

## I Introduction

The quasi-one-dimensional (quasi-1D) $`S=1/2`$ antiferromagnet BaCu2Si2O7 is recognized as an almost ideal model material for studying exotic spin dynamics in weakly-interacting quantum spin chains.Tsukada et al. (1999) In a series of papersZheludev et al. (2000); Kenzelmann et al. (2001); Zheludev et al. (2002) we have demonstrated that in this compound it is possible to observe a separation between single-particle spin wave excitations, characteristic of a classical system with long-range magnetic order, and a gapped continuum of states that is a property of the quantum 1D Heisenberg model. Quite recently however, another aspect of the physics of BaCu2Si2O7 attracted a great deal of attention. It was found that in the ordered state of BaCu2Si2O7, an application of a magnetic field along the direction of ordered moment leads to two consecutive spin re-orientation transitions.Tsukada et al. (2001) This phenomenon is in contrast with the typical behavior of a conventional antiferromagnet in a field, where only a single “spin-flop” transition occurs.

By far, the largest energy scale of magnetic interactions in BaCu2Si2O7 is that of antiferromagnetic nearest-neighbor coupling within the spin chains, that run along the $`c`$ axis of the orthorhombic crystal structure (space group Pnma, $`a=6.862(2)`$Å, $`b=13.178(1)`$Å and $`c=6.897(1)`$Å.Oli ). The in-chain exchange constant is $`J=24.1`$ meV, as previously determined from inelastic neutron scattering experiments and bulk susceptibility measurements.Tsukada et al. (1999); Kenzelmann et al. (2001) Exchange coupling between the chains is much weaker, with $`J_x=0.46`$ meV and $`J_y=0.20`$ meV, for nearest-neighbors along the $`a`$ and $`b`$ axes, respectively, and $`J_3=0.15`$ meV along the $`(110)`$ diagonal.Kenzelmann et al. (2001) The spin wave spectrum shows a small anisotropy gap of about 0.3 meV, corresponding to an easy axis along the $`(001)`$ crystallographic direction.Kenzelmann et al. (2001) Long-range magnetic ordering occurs at a low Néel temperature of $`T_N9.2`$ K.Tsukada et al. (1999) The ordered magnetic moment is rather small, $`m_0=`$0.15 $`\mu _\mathrm{B}`$ per Cu2+.Zheludev et al. (2000); Kenzelmann et al. (2001) In the ordered state the spins are parallel to the magnetic easy axis. Nearest-neighbor spins along the $`b`$ and $`c`$ axes are aligned antiparallel relative to each other, while nearest neighbors along $`a`$ are ferromagnetically aligned.Tsukada et al. (1999); Kenzelmann et al. (2001) The spin structure is visualized in Fig. 1. The phase transitions that are the subject of this work occur at relatively low fields $`H_{c1}=2.0`$ T and $`H_{c2}=4.9`$ T, applied along the easy axis, and produce characteristic jumps in the magnetization curve.Tsukada et al. (2001) To date, the spin structure in the high-field phases has not been determined. However, in Ref. Tsukada et al., 2001 it was suggested that upon going through the first transition, the spins “flop” into the $`(a,b)`$ crystallographic plane. The second transition was then described as an unusual spin rotation within the plane perpendicular to the field, and is somewhat similar to that recently found in K2V3O8.Lumsden et al. (2001)

Based on an analysis of the high-field magnetization data, it was conjecturedTsukada et al. (2001) that the unusual second spin re-orientation in BaCu2Si2O7 is a result of competition between Dzyaloshinskii-Moriya (DM) off-diagonal exchange interaction in the spin chains and isotropic exchange coupling between the chains. In the present work we report the results of magnetic neutron diffraction studies of the various spin-flop states in BaCu2Si2O7. We find the results to be inconsistent with the previously proposed model. It appears that inter-chain coupling dominates over DM interactions, and that the series of spin re-orientation transitions can not result from a competition between the two effects alone. Alternative mechanisms for the two-stage spin flop transition in BaCu2Si2O7 are discussed.

## II Experimental

A $`5\times 5\times 4`$ mm3 sample was mounted in a 6 T cryomagnet on the D23 lifting counter diffractometer at Institut Laue Langevin. The $`(0,0,1)`$ axis was aligned parallel to the field direction. With 8 magnetic Cu2+ per unit cell, in the ordered state all magnetic reflections coincide with nuclear ones. In addition, the ordered magnetic moment in BaCu2Si2O7 is small, and detecting the magnetic contributions to Bragg intensities is rather difficult. The measurements were therefore performed in differential mode. At $`H=0`$ T, $`4`$ T (below $`H_{c2}`$) and 5 T (above $`H_{c2}`$) the diffractometer was consecutively positioned on each accessible Bragg reflection. Without moving any instrument motors, the peak intensity was measured first at 2 K and then again at 10 K. This measurement yielded the ratio of nuclear and magnetic scattering intensities. The absolute values of nuclear intensities were separately measured at 10 K in standard rocking scans. The rocking curves also provided an estimate of the background around each Bragg position. The measured integrated nuclear intensities, backgrounds, and peak intensity ratios were then used to estimate integrated magnetic intensities. In another experiment, the peak intensities of $`(\overline{4}1\overline{1})`$, $`(0\overline{3}\overline{1})`$ and $`(0\overline{5}\overline{1})`$ reflections were monitored as a function of applied field in the range 0–5.5 T. Peak intensities measured in this fashion were assumed to be proportional to integrated intensities, with the normalization factor determined at $`H=0`$, as described above.

## III Results

Due to the strong in-chain coupling, any canting of neighboring spins away from a perfect antiferromagnetic alignment within each chain will be small and undetectable in our measurements. In analyzing the data, we could therefore safely rely on an approximation in which within each chain all spins are collinear. Inter-chain coupling is considerably weaker than the in-chain one, and is comparable in strength to the Zeeman energy associated with experimental fields. For this reason the relative orientation of staggered moments in neighboring chains could not be considered fixed, as it may actually vary with field. Since there are four chains within each unit cell of BaCu2Si2O7, solving the magnetic structure at each field was reduced to determining the orientations of four spins, one from each chain, labeled as $`𝐒_1`$$`𝐒_4`$ in Fig. 1.

### III.1 Spin structure at zero applied field: Phase I

As a first step, we have verified that our new data collected at $`H=0`$ are consistent with the structure previously reported in Ref. Kenzelmann et al., 2001. The spins are parallel to the $`c`$-axis, with nearest neighbor spins along the $`c`$ and $`b`$ axes antiparallel to each other, and a parallel alignment of nearest neighbor spins along the $`a`$ axis, as shown in Fig. 1. This type of spin arrangement minimizes all inter-chain interaction energies that were independently measured using inelastic neutron scattering.Kenzelmann et al. (2001) Magnetic Bragg intensities calculated for this $`c`$-axis colinear state are in good agreement with those measured at $`H=0`$ in the present work. This comparison is made in Table 1.

### III.2 Spin structure at $`H=4`$ T: Phase II

The magnetic structure at $`H=4`$ T was solved using a reverse Monte Carlo analysis, which unambiguously pointed to a planar state with all spins in the $`(a,b)`$ plane and a canting by an angle $`\varphi 80^{}`$, as shown in Fig. 2 on the left. This structure is a slight distortion of a colinear one, with all spins parallel to the the $`b`$ axis, and the same relative alignment of spins in adjacent chains as at $`H=0`$. A least-squares refinement of the canting angle and the ordered moment then yielded $`\varphi =78(3)^{}`$, $`m_0=0.18(2)`$ $`\mu _\mathrm{B}`$, with a residual $`\chi ^2=2.4`$.

A comparison between measured and calculated magnetic intensities is given in Table 2. Assuming $`\varphi =90^{}`$ ($`b`$-axis colinear state), one can still get a reasonably good agreement with experiment, with $`\chi ^2=2.8`$. Note that the magnetic structure of Phase II proposed in Ref. Tsukada et al., 1999 (all spins parallel or antiparallel to the $`a`$-axis) is totally inconsistent with the present experiment, corresponding to $`\chi ^2=64`$, even when $`m_0`$ is optimized to best fit the data.

### III.3 Spin structure at $`H=5`$ T: Phase III.

Fewer data points were collected at $`H=5`$ T than for the other two field values (see Table 3). The best agreement with the data was obtained assuming a collinear structure with all spins pointing along the $`a`$ axis (Fig. 2, right), and the same relative alignment of spins in adjacent chains as at $`H=0`$. The refined value for the ordered moment is $`m_0=0.17(2)`$ $`\mu _\mathrm{B}`$, and for this model $`\chi ^2=4.5`$. Again we note that the spin arrangement proposed in Ref. Tsukada et al., 2001 for the high-field phase (all spins along the $`b`$ axis and nearest-neighbor spins in the $`a`$ direction aligned parallel to each other) is totally inconsistent with experiment, corresponding to $`\chi ^2=158`$.

## IV Field dependence.

The measured field dependencies of $`(\overline{4}1\overline{1})`$, $`(0\overline{3}\overline{1})`$ and $`(0\overline{5}\overline{1})`$ reflections are shown in Fig. 3. The data for Phases II and III were analyzed using the established models. The value of magnetic moment $`m_0`$, as well as the canting angle $`\varphi `$ for Phase II were refined to best-fit the three intensities at each field. The results of the fit are shown in solid lines in Fig. 3, and the field dependence of the fit parameters is plotted in Fig. 4.

To obtain a good fit to the field dependencies measured in Phase I, an additional parameter had to be introduced to account for the large change in intensities just below the phase transition at $`H_{c1}`$. By trial and error, it was found that the best way to reproduce this effect is to assume that the spin structure, as a whole, tilts towards the $`b`$ axis (Fig. 5). In this model, when the direction of staggered moment passes through the $`(0\overline{3}\overline{1})`$ direction, the corresponding Bragg intensity shows a characteristic dip (arrow in Fig. 3a). This behavior is a result of the neutron polarization intensity factor going to zero for scattering vectors that are parallel to ordered magnetic moment. A smaller dip, at a slightly larger field is seen in the $`(0\overline{5}\overline{1})`$ peak as well (arrow in Fig. 3b). By varying $`m_0`$ and the tilt angle $`\alpha `$, a good fit to the data collected at $`0<H<H_{c1}`$ could be obtained. The result of the fit is shown in solid lines in Fig. 3, and the field dependence of the fit parameters is plotted in Fig. 4.

Note that the angle $`\alpha `$ is shown as negative at small fields. In fact, the tilt $`\alpha `$ probably remains close to zero and positive, except immediately before the transition, where it rapidly increases to $`90^{}`$. The negative values obtained in the fit are an artifact of using only three measured intensities in the analysis. For small $`\alpha `$ the experiment is fairly insensitive to this parameter, due to a $`\mathrm{cos}^2\alpha `$ term in the intrinsic polarization factor for the scattering cross section and the fact that all data were collected close to the $`(001)`$ reciprocal-space plane. However, as $`\alpha `$ approaches $`90^{}`$, the polarization sensitivity is greatly enhanced.

## V Discussion

The data analysis described above was based on the known room-temperature crystal structure. The use of room-temperature structure factors to determine spin orientations at low temperatures may lead to systematic errors. Additional errors may be associated with measuring magnetic intensities in the temperature-differential mode. Nevertheless, the results allow us to establish the main features of the three phases: 1) Nearest-neighbor spins within each chain remain almost anti-parallel to each other at all times. 2) The relative alignment of adjacent spins from neighboring chains does not deviate significantly from that at $`H=0`$, with ferromagnetic spin alignment along $`a`$ and antiferromagnetic along $`b`$. Even in Phase II the canting is not very dramatic and vanishes just before the transition at $`H_{c2}`$. 3) The spins are oriented roughly along the $`c`$, $`b`$ and $`a`$ axis in phases I, II and III, respectively. 4) The ordered moment is only slightly field-dependent and remains close to $`0.15\mu _\mathrm{B}`$.

In weakly coupled quantum spin chains the ordered moment $`m_0`$ is proportional to $`\sqrt{|J_{}|/J}`$, where $`J_{}`$ is the effective mean-field inter-chain coupling strength.Sch  At $`H=0`$ the structure is such that $`|J_{}|=J_b+2J_3+|J_a|`$.Kenzelmann et al. (2001) The fact that $`m_0`$ remains practically constant when going through the phase transitions, means that in Phases II and III, $`J_{}`$ does not change much and nearest-neighbor spins remain close to being parallel to one another along the $`a`$ axis and antiparallel along the $`b`$ and $`c`$ directions. This can be seen as an important self-consistency check for our determination of the spin orientations in the three phases. The overall gradual increase of the ordered moment as a function of $`H`$ seen in Fig. 4a is to be expected, and is due to the suppression of 1D quantum spin fluctuations by the symmetry-breaking external field. The slight suppression of $`m_0`$ observed in Phase II (Fig. 4a) is a result of the canted geometry, in which $`|J_{}|`$ is reduced compared to that in Phase I. In the presence of canting ($`\varphi <90^{}`$) the effective inter-chain coupling is given by $`|J_{}|=2J_3(J_b+|J_a|)\mathrm{cos}(2\varphi )`$. In particular, just above $`H_{c1}`$, $`\varphi =60^{}`$, and the drop in $`m_0`$ expected on going through the first phase transition due to the abrupt change in $`|J_{}|`$ is about 23%. This value is in reasonable agreement with the observed $`15\%`$ effect.

Given the available data, the spin re-orientation at $`H_{c1}`$ can be described as a more or less conventional spin-flop transition, that one expects for a classical antiferromagnet in a magnetic field applied along the easy axis. In zero field, the easy-$`c`$-axis anisotropy energy is minimized when all spins are perpendicular to the $`(a,b)`$ plane. To take advantage of the Zeeman energy in applied fields, the spins have to align themselves perpendicular to $`c`$, to allow a canting of the rigid antiferromagnetic structure in the direction of $`𝐇`$. Dzyaloshinskii-Moriya interactions and inter-chain coupling may modify the transition field and produce the small canting found in Phase II, but do not affect the basic physics of the effect. The fact that the transition appears to have a precursor, the spin structure tilting by an angle $`\alpha `$ at $`H<H_{c1}`$, may have several explanations. One possibility is that this phenomenon is result of an imperfect alignment of the external magnetic field relative to the $`c`$ axis of the crystal. Alternatively, it can be an intrinsic effect, due to the presence of Dzyaloshinskii-Moriya interaction, or to the existing structural canting of the local anisotropy axes of Cu2+ relative to the $`c`$ axis. To fully resolve this problem, a careful determination of the magnetic structure just below and just above $`H_{c1}`$ will be required.

The second phase transition is much more unusual, since it involves a spin rotation in the plane perpendicular to $`𝐇`$. It is clear that such behavior requires the presence of off-diagonal exchange interactions in the system. As recently observed in K2V3O8, for example, a spin rotation around the field direction can be caused by a competition between DM interactions and magnetic anisotropy.Lumsden et al. (2001) This type of behavior is known since earlier studies of hematite.hem  For BaCu2Si2O7 though, a totally different explanation was proposed in Ref. Tsukada et al., 2001, and the transition was attributed to a competition between DM coupling in the chains and isotropic inter-chain interactions, a mechanism similar to that previously discussed in relevance to the anomalous spin-flop behavior of La2CuO4.Thio et al. (1990)

The key argument of Ref. Tsukada et al., 2001 is that the main components of the Dzyaloshinskii vector for the nearest-neighbor bond within the chains lie in the $`(a,b)`$ plane, and alternate sign from one bond to the next. In this case, DM interactions produce a weak-ferromagnetic canting of the spins within each chain even in zero applied field. In the presence of an external field, additional canting is due to Zeeman energy. The free energy of the system is minimized when the two canting effects are in the same direction. For this to be the case, nearest-neighbor spins along the $`a`$ axis have to be almost antiparallel to each other, since Dzyaloshinskii vectors in the adjacent chains are antiparallel as well. This preferred antiparallel orientation is in competition with ferromagnetic inter-chain coupling along the $`a`$ axis. In the model proposed in Ref. Tsukada et al., 2001, the frustration is resolved in Phase II by having all spins aligned almost parallel to the Dzyaloshinskii vector (roughly along the $`a`$ axis), which practically eliminates the Dzyaloshinskii energy altogether, and the ferromagnetic inter-chain coupling energy along the $`a`$ axis is minimized. In Phase III the model predicts that the combination of DM and Zeeman energies wins over inter-chain interactions: the spins rotate to be almost parallel to the $`b`$ axis and cant in the direction of the applied field to minimize both energies. At the same time, nearest-neighbors along the $`a`$ axis become almost antiparallel, despite a ferromagnetic coupling in this direction.

The results presented above clearly demonstrate that the actual spin arrangements in BaCu2Si2O7 in Phases II and III are totally different from those previously conjectured. In Ref. Tsukada et al., 2001 it was emphasized that the Dzyaloshinskii vectors not being strictly parallel to the $`a`$ axis, and additional inter-chain interactions along the $`(0,1,0)`$ and $`(1,1,0)`$ directions being present, the critical field and the details of the spin structures may deviate from those predicted by the simplified model. However, even in the more general case, if the $`H_{c2}`$ transition in BaCu2Si2O7 was indeed driven by a competition between DM interactions and isotropic inter-chain coupling, the energy of the latter would have to change abruptly upon going through the transition. The present data clearly show that this is not the case, since at $`H_{c2}`$ the spin structure rotates as a whole, and the relative spin orientations, and with them the inter-chain exchange energies, remain practically unchanged. It appears that inter-chain coupling dominates over DM interactions in all three phases. Because of the canting in Phase II, the energy of inter-chain interactions is actually better minimized in Phase III, rather than Phase II. This sequence is opposite to that emerging from the model proposed in Ref. Tsukada et al., 2001.

We believe that at the present stage it is unwise to attempt a quantitative explanation of the phase transitions in BaCu2Si2O7 based on some specific spin Hamiltonian. Too many parameters are involved to make such an analysis unambiguous. The energy scale of inter-chain interactions and $`c`$-axis anisotropy are indeed similar to that defined by transition fields. However, the Hamiltonian is expected to have other terms on the same energy scale, and these can not be dismissed. In particular, measurements of the spin wave spectrum contain evidence of an in-plane anisotropy of the magnitude $`0.15`$ meV. This anisotropy is most likely due to the so-called Kaplan Schekhtman Entin-Wohlman Aharony (KSEA) interactions.KSE  The KSEA term is a companion to DM interactions, and is an easy axis parallel to the Dzyaloshinskii vector. KSEA interactions are believed to be the driving force in the spin reorientation transition in K2V3O8.Lumsden et al. (2001) The similarity of the behavior of the latter compound with that found in BaCu2Si2O7 may indicate that anisotropy effects, and KSEA interactions in particular, could be responsible for the $`H_{c2}`$ transition in the silicate as well.

## VI Conclusion

In summary, the actual mechanism of the exotic two-stage spin flop transition in BaCu2Si2O7 remains a mystery. What is clear though, is that DM interactions and inter-chain exchange coupling alone can not sufficient to explain phenomenon. It appears more likely that both transitions in BaCu2Si2O7 are caused by a competition between DM interactions and diagonal magnetic anisotropy effects, of which more needs to be learned.

###### Acknowledgements.

One of the authors (I. T.) thanks J. Takeya for useful discussions. This work is supported in part by the Grant-in-Aid for COE Research “SCP coupled system” of the Japanes Ministry of Education, Science, Sports, and Culture. Oak Ridge National Laboratory is managed by UT-Battelle, LLC for the U.S. Department of Energy under contract DE-AC05-00OR22725.


text:
# B1422+231:The influence of mass substructure on strong lensing

## 1 Introduction

Gravitational lens systems with multiply imaged quasars are an excellent tool for studying the properties of distant galaxies. In particular, they provide the most accurate mass measures for the lensing galaxy. Besides the mass profiles, one can also gain information about evolution (Kochanek et al. 2000) and extinction laws (Fassnacht et al. 1999). Strong lensing is also a very promising and robust tool to measure the Hubble constant (Refsdal 1964). The success of this method, however, depends strongly on how well the mass model is constrained.

It turns out that image positions can be fit quite accurately with simple, smooth elliptical models (Keeton et al. 1997). Since the number of observational constraints from image positions is small, one wants to include the flux information. The optical fluxes, however, should not be used, as they might be affected by microlensing and/or dust obscuration (Chang & Refsdal 1979). Radio fluxes, on the other hand, can provide further constraints in lens modelling.

Fitting the fluxes very accurately turns out to be difficult in many lens systems. Models for MG0414+0534 (Falco et al. 1997; Keeton et al. 1997), PG1115+080 (Keeton & Kochanek 1997) and B1422+231 (e.g. Kormann et al. 1994) all show the same failure, namely the observed flux ratios are very different from what one would expect for the image configurations from a smooth model. In particular, the gravitational lens system B1422+231 was explored in detail, and as was first mentioned by Mao & Schneider (1998), mass substructure in the lens galaxy might provide an explanation for the failures in flux modelling.

A question arises as to whether there is something special about these quadruple systems, or does the discrepancy simply arise due to the fact that our smooth models are oversimplified? In other words, we are asking how “well” (for the purpose of strong lensing) can the smooth models used to fit the data represent a real galaxy? N-body simulation data can provide a realistic description of a galaxy mass distribution, thus giving a possibility to probe its effect on strong lensing.

In this paper, we will study the influence of the substructure on the lens system B1422+231, using VLBI radio measurements of the system by Patnaik et al. (1999). In Sec. 2 we give a description of the lens system and data used. The method is outlined in Sect. 3 and the results on fitting the system with smooth mass model are presented. In Sect. 5 the model accounting for the substructure is presented and in Sect. 5.2 the deconvolved image shape information is added to the fit. Sect. 6 gives the description of the method used to investigate lensing by an N-body simulated galaxy. We describe how we obtained synthetic data of four image systems, and the results of fitting such systems with lens models are presented. Finally we draw some conclusions in Sect. 7.

This work is an abbreviated version of Bradač (2001). In the course of writing this paper, several related papers on substructure of lens galaxies have been submitted (Metcalf & Madau 2001; Chiba 2001; Dalal & Kochanek 2001; Metcalf & Zhao 2001; Keeton 2001). In Metcalf & Madau (2001) and Chiba (2001), the authors use semi-analytical descriptions to account for mass clumps typical for N-body simulations in a statistical fashion. In addition, Chiba (2001) tests his prediction on two known systems with four images. Dalal & Kochanek (2001), Metcalf & Zhao (2001), and Keeton (2001) predict the properties of substructure needed to constrain lens systems and compare their results with the CDM predictions. In addition, Keeton (2001) accounts for the difference between radio and optical fluxes by investigating radio and optical sources of different sizes. The present work is differs from the afore mentioned in that slightly different models for the lens in B1422+231 system are used and that we also include image shape constraints in the fit. Further, we investigate four image systems generated by using a CDM N-body simulated galaxy, rather than an analytic approximation for the statistics of mass-substructure.

## 2 The mystery of B1422+231

The gravitational lens system B1422+231 was discovered in the course of the JVAS survey (Jordell Bank – VLA Astrometic Survey) by Patnaik et al. (1992). It consists of four image components. The three brightest images A, B, and C (as designated by Patnaik et al. 1992) are fairly collinear. The radio flux ratio between images A and B is approximately 0.9, while image C is fainter (flux ratio C to B is approximately 0.5). Image D is further away and is much fainter than the other images (with flux ratio D:B of 0.03). The most recent available radio data for the image positions and fluxes were obtained from the polarisation observations made at $`8.4\mathrm{GHz}`$ using the VLBA and the 100m telescope at Effelsberg from Patnaik et al. (1999) and are listed in Table 1. For each of the components, the authors measured positions (relative to the image B) and fluxes as well as the deconvolved image shapes. Here and through the paper we are using a notation where $`(\theta _1,\theta _2)`$ are the angular coordinates in the lens plane and $`(\beta _1,\beta _2)`$ in the source plane. $`\theta _1`$ and $`\beta _1`$ increase in the negative RA direction.

The radio source of this lens system is associated with a $`15.5\mathrm{mag}`$ quasar at a redshift of $`3.62`$ (Patnaik et al. 1992). The lensing galaxy has been observed in the optical; its redshift has been determined to be $`0.338`$ and its position relative to image B has been measured (Impey et al. 1996). The main lens galaxy is a member of a compact group with a median projected radius of $`35h^1\mathrm{kpc}`$ and velocity dispersion of $`550\mathrm{km}\mathrm{s}^1`$ (Kundic et al. 1997).

Several groups have tried to model B1422+231 (Hogg & Blandford 1994; Kormann et al. 1994; Keeton et al. 1997; Mao & Schneider 1998) and all of them have experienced difficulties in fitting it. As we used data with even more precise image positions one might expect that it would become even harder to model the system. However, as already pointed out by some authors, the difficulties do not lie in fitting the image positions but rather in the flux ratios.

It turns out that simple, smooth models fail to reproduce the radio as well as the optical flux ratios of the system. While the mismatch in optical data might still be due to the microlensing and/or dust obscurations, this can probably not explain why such models are not successful when fitting radio flux ratios. Mao & Schneider (1998) have proposed a lens model that accounts for the substructure in the lensing galaxy. They concluded that one needs a surface mass density perturbation of the order of 1 % of the critical surface mass density in order to change the flux ratios to the observed values.

In order to succeed in fitting the image flux ratios, one needs to consider more sophisticated models. However, such models also require the use of additional parameters. Therefore it is very difficult to ensure a constrained model that accounts for the substructure using as constraints only image positions, flux ratios, and the galaxy position. For this reason we also included the axis ratios and the orientation angles of the deconvolved images as additional constraints. They were obtained from Patnaik et al. (1999) and are given in Table 2. Listed are the absolute values of ellipticities and the orientation angles of the fitted elliptical Gaussians together with uncertainties.

From the configuration and fluxes of the four images we can gain some qualitative constraints on the lens model. Image D is much fainter than the rest which can mean that it is either highly demagnified or that the other three images are highly magnified. Because the position of the primary lens galaxy is known, the possibility of image D being highly demagnified is ruled out. Namely, image D does not lie at a position of high surface mass density and thus high demagnification. Images A, B, and C are therefore highly magnified. In order to get three highly magnified images with a fairly collinear configuration, the source has to lie close to (and inside) a cusp (e.g. Schneider et al. 1992 and references therein).

For a source position sufficiently close to and inside a cusp there exists a relation which states that the sum of the fluxes of the outer two magnified images (in our case images A and C) is the same as the flux of the middle image – B (e.g. Mao 1992; Schneider & Weiss 1992). This rule is strongly violated in the lens system B1422+231. Actually, one can get a deviation from this relation if the source is away from the cusp or if there is substructure in the system (i.e. fluctuations on a scale smaller than the separations between A, B, and C). Also the directions of elongation of the images indicate that the source is indeed close to a cusp and thus argue in favour of substructure in the system (Mao & Schneider 1998).

## 3 Lens modelling

First we considered two standard gravitational lens models since their application to the Patnaik et al. (1999) data has not yet been discussed in the literature.

The standard approach to model a strong lens system is to define the goodness-of-fit function $`\chi ^2`$, a measure of the deviation of the predicted and observed image properties. We perform the minimisation in the image plane. The $`\chi ^2`$ for the positions of the images reads

$$\chi _{\mathrm{pos}}^2=\underset{i=1,3,4}{}\frac{\left|\mathrm{\Delta }𝜽^{i,\mathrm{obs}}\mathrm{\Delta }𝜽^{i,\mathrm{mod}}\right|^2}{\sigma _{i,\mathrm{pos}}^2},$$
(1)
where $`\mathrm{\Delta }𝜽^{i,\mathrm{obs}}`$ and $`\mathrm{\Delta }𝜽^{i,\mathrm{mod}}`$ are the observed and the modelled position of the $`i`$-th image relative to a chosen origin (in our case image B, denoted by $`i=2`$), respectively. Note that the sum extends only over the images A, C, and D.

To be more precise, $`\mathrm{\Delta }𝜽^{i,\mathrm{mod}}`$ is obtained as

$$\mathrm{\Delta }𝜽^{i,\mathrm{mod}}=𝜽^{i,\mathrm{mod}}𝜽^{2,\mathrm{mod}},$$

where the vectors $`𝜽^{i,\mathrm{mod}}`$ and $`𝜽^{2,\mathrm{mod}}`$ are measured with respect to an independent reference point. The same is true for $`\mathrm{\Delta }𝜽^{i,\mathrm{obs}}`$,

$$\mathrm{\Delta }𝜽^{i,\mathrm{obs}}=𝜽^{i,\mathrm{obs}}𝜽^{2,\mathrm{obs}},$$

which is the observed image position as listed in Table 2. Therefore when the image positions of A, C, and D w.r.t. the image B are measured, they all contain the uncertainty of the origin (image B). Thus, the uncertainties of the measurements of each of the image positions are correlated and we are forced to calculate the $`\chi _{\mathrm{pos}}^2`$ as in (1). We therefore denote $`\sigma _{i,\mathrm{pos}}`$ to be the uncertainty of the measured relative image position.

Similarly, we define for the galaxy position

$$\chi _{\mathrm{galaxy}}^2=\frac{\left|\mathrm{\Delta }𝜽_{\mathrm{lens}}^{\mathrm{obs}}\mathrm{\Delta }𝜽_{\mathrm{lens}}^{\mathrm{mod}}\right|^2}{\sigma _{\mathrm{gal},\mathrm{pos}}^2},$$
(2)
with $`\sigma _{\mathrm{gal},\mathrm{pos}}`$ being the uncertainty of the relative lens position. The contribution to the $`\chi ^2`$-function from the fluxes is given by

$$\chi _{\mathrm{flux}}^2=\underset{i=1}{\overset{4}{}}\frac{\left(F^{i,\mathrm{obs}}F^{i,\mathrm{mod}}\right)^2}{\sigma _{i,\mathrm{flux}}^2},$$
(3)
where $`F^{i,\mathrm{mod}}`$ denotes the flux of the $`i`$-th image obtained from the model and $`F^{i,\mathrm{obs}}`$ is the observed flux.

Further we can use deconvolved image shapes as constraints of the model. We work in terms of the complex ellipticity,

$$ϵ_i:=\left|ϵ_i\right|\mathrm{e}^{2\mathrm{i}\phi _i}.$$

Each image is thus described by the absolute value of ellipticity $`\left|ϵ_i\right|`$ and the corresponding position angle $`\phi _i`$ (measured w.r.t. $`\theta _1`$-axis). From the given data the absolute value of ellipticity is calculated as

$$\left|ϵ_i\right|=\frac{a_ib_i}{a_i+b_i},$$
(4)
where $`a_i`$ and $`b_i`$ are major and minor semi axes of the fitted elliptical Gaussian. Two additional terms,

$$\chi _ϵ^2=\underset{i=1}{\overset{4}{}}\frac{\left(\left|ϵ_i^{\mathrm{obs}}\right|\left|ϵ_i^{\mathrm{mod}}\right|\right)^2}{\sigma _{\left|ϵ\right|,i}^2},$$
(5)
and

$$\chi _\phi ^2=\underset{i=1}{\overset{4}{}}\frac{\left(\phi _i^{\mathrm{obs}}\phi _i^{\mathrm{mod}}\right)^2}{\sigma _{\phi ,i}^2}.$$
(6)
have to be added to the $`\chi ^2`$-function. For simplicity we assumed the errors on ellipticity and position angle to be uncorrelated. The $`\chi ^2`$ function we want to minimise is simply given by the sum

$$\chi _{\mathrm{tot}}^2=\chi _{\mathrm{pos}}^2+\chi _{\mathrm{galaxy}}^2+\chi _{\mathrm{flux}}^2+\chi _ϵ^2+\chi _\phi ^2.$$
(7)
In order to obtain the values of $`\mathrm{\Delta }𝜽^{i,\mathrm{mod}}`$ from the lens equation (see e.g. Schneider et al. 1992), we used the Numerical Recipes MNEWT routine from Press et al. (1992) for solving a set of two non-linear equations. The $`\chi ^2`$-function was then minimised with respect to the model parameters using POWELL, a multi-dimensional minimisation routine, also from Press et al. (1992).

One of the consideration one faces when doing lens modelling is the uniqueness of the solution. Minimisation methods such as POWELL do not give an answer as to whether a minimum found is a local or global one. One can partly solve this by using simulated annealing routines, such as AMOEBSA from Press et al. (1992). However, these routines are much more CPU-time consuming, and they are also not completely foolproof in finding an absolute minimum.

As an additional check, to minimise the chance of obtaining a secondary minimum as a result, we have run POWELL for several different initial conditions. Unless the global minimum is extremely narrow in parameter space (which is rather unlikely) we can be confident of obtaining the global minimum.

## 4 Modelling B1422+231 with smooth models

We used a singular isothermal ellipsoid with external shear from Kormann et al. (1994) (hereafter SIE+SH) and a non-singular isothermal ellipsoid model with external shear (NIE+SH) from Keeton & Kochanek (1998) to fit the image positions and fluxes of B1422+231. The explanation of the model parameters for the models we used are given in Table 3.

We have applied the fitting procedure described above to the radio data, using image positions, fluxes, and their uncertainties from Patnaik et al. (1999), listed in Table 1. The optical position of the galaxy was taken from Impey et al. (1996). Although the image positions are very accurate (of the order of $`50\mu \mathrm{arcsec}`$), we have no difficulties fitting them, and so the $`\chi ^2`$ contribution from the image positions drops to zero (see Table 4). However, as already pointed out in previous works on B1422+231, the model completely fails in predicting the image fluxes. In particular image A is predicted too dim (the modelled flux ratio A:B turned out to be 0.80, much below the measured value of 0.93). We have also tried to model the system with a NIE+SH model; however, the $`\chi ^2`$ did not improve significantly.

Next we disregarded the flux of image A as a constraint, trying to see whether we can get a good fit for the rest of the images. Reducing the number of constraints by one, we are left with one degree of freedom. The results of the fitting are given in Table 4. The fit is still not satisfactory, thus indicating that changes of magnifications are needed not only for image A.

The choice of these two particular macromodels was due to their simplicity and analytic solutions for the deflection angle. Of course these are not the only models one can use, and one might hope to explain the discrepant flux ratios by using another family of models. In order to effectively change the flux ratios, one has to change the magnification gradient. The simplest way of changing a gradient is to change the power-law behaviour of the potential.

Although the success of such modelling is very unlikely (as shown by Mao & Schneider 1998), we have tried to use the SIE+SH model, adding an additional black hole with an Einstein radius fixed to $`10^2`$ of the main galaxy Einstein radius, centred on the lens galaxy position. The corresponding mass of the black hole is $`10^8\mathrm{M}_{}`$, in agreement with the correlation of the black hole mass with the velocity dispersion of its host bulge (Ferrarese & Merritt 2000). Such a model is plausible, however we found no significant improvement of the fit (see Table 4 for results). Also, considering the black hole mass as a parameter did not significantly improve the outcome.

Although other smooth models can still be investigated, it seems unlikely that another smooth model can explain all four image fluxes simultaneously. Even when one disregards the flux of image A, a smooth macromodel seems to be incapable of explaining the remaining flux ratios.

## 5 Models with substructure

In the previous section it turned out that A:B flux ratio causes the biggest difficulty in fitting the B1422+231. Since the radio and optical flux ratios are very different, one is tempted to exclude it from the $`\chi ^2`$ measure (see Mao & Schneider 1998, Chiba 2001).

However, one can also try to deal with this problem in another way. Adding a small perturber at the same angular diameter distance as the primary lens and at approximately the same position as image A can change the flux ratio A:B substantially. On the other hand, calculations show (Mao & Schneider 1998) that such a perturber does not affect the positions of any of the images appreciably. Furthermore, a small perturber can also change the flux ratio of the other images slightly and this might help to improve the results from the previous section.

We model the perturber as a non-singular isothermal sphere. This gives two additional parameters (two perturber positions, since we keep the line-of-sight velocity dispersion and core radius of the perturber fixed) to the macro (SIE+SH) model. The choice of the perturber being modelled as non-singular is due to the fact that a singular isothermal sphere (with the same Einstein radius) is more likely to give rise to additional (observable) images.

We are aware of the fact that the choice of this particular model for the substructure is oversimplified in many ways. However, we are not trying to constrain the nature of substructure in this case; which is impossible due to the number of constraints available. In the following we investigate whether the peculiarity of the fluxes can be explained by a small satellite galaxy in the neighbourhood of one (or two) of the images.

### 5.1 Modelling B1422+231 using the substructure model

For modelling B1422+231 we fixed the line-of-sight velocity dispersion of the NIS perturber to $`\sigma _{\mathrm{NIS}}=10\mathrm{km}\mathrm{s}^1`$ (equivalent to an Einstein radius of approximately $`2\mathrm{mas}`$) and its core radius to $`\theta _\mathrm{c}=20\mathrm{mas}`$. A perturber with these properties does not affect the image positions significantly; on the other hand it can substantially change the magnification at the position of one of the images.

When fixing the core radius one has to be aware that not only an SIS, but also an NIS lens with small enough core radius might give additional images. Therefore, an NIS perturber should have a core radius much bigger than the Einstein radius, in order not to produce additional images. We have checked that indeed no additional observable images are predicted by the model.

For this purpose we define the function $`f`$ on a grid of points $`𝜽_j`$

$$f=\left|𝜷(𝜽_j)𝜷_\mathrm{s}\right|^2,$$
(8)
where $`𝜷(𝜽_j)`$ is the calculated source position corresponding to the point $`𝜽_j`$. The position $`𝜷_\mathrm{s}`$ is defined as the average position for the source predicted by the four observed image positions

$$𝜷_\mathrm{s}=\frac{1}{4}\underset{i=1}{\overset{4}{}}𝜷\left(𝜽^{i,\mathrm{obs}}\right).$$
(9)
The function $`f`$ vanishes only around the positions where images are observed, provided that the grid resolution is high enough. In our case the grid resolution was $`0.05\mathrm{mas}`$ (small compared to the Einstein radius of the perturbing galaxy, $`\theta _{\mathrm{E},\mathrm{NIS}}=3\mathrm{mas}`$). We indeed found only four such regions, and they correspond to the four observed images.

The resulting model has 12 parameters, which leaves us 0 degrees of freedom. The $`\chi ^2`$ has decreased by a factor of more than 20 compared with the SIE+SH model; however, since we have zero degrees of freedom we expect $`\chi ^2`$ to vanish if the model is realistic (and if the $`\chi ^2`$ technique is an adequate method). The family of models considered thus does not seem to be adequate for the description of the galaxy in B1422+231 lens system.

We again see that the flux ratio of images A:B is not the only problem when dealing with fluxes. Adding a small perturber close to image A allows us to adjust the modelled flux of A such, that it does not give any contribution to the $`\chi ^2`$-function. Still the remaining image fluxes are not fit perfectly, leading to a possible conclusion that all images are affected by the mass substructure.

### 5.2 Using image shapes as constraints

In order to increase the number of degrees of freedom we are dealing with, further constraints have to be added. We therefore used the deconvolved image shapes (see Table 2). Due to the problems in determining the image shapes we describe later, we will use them as constraints only in this section.

The $`\chi ^2`$ minimisation was done according to the previous section, and the results are presented in Table 5.

We have decreased the value of the core radius of the perturber in order to get higher magnification gradients, since large changes in image shapes as predicted by a smooth model was needed. Again, no additional observable images are produced by the perturber(s).

In the Patnaik et al. (1999) paper the uncertainties on the image shapes are not listed. The image shapes are obtained by fitting Gaussian profiles to the map, and then deconvolved using the known beam-shape. The uncertainties are therefore just a rough estimate, since one can not quantitatively account for the error of such fitting. In fact, all the images exhibit non-Gaussian features a Gaussian model is an oversimplification for the image shape description.

The resulting $`\chi ^2`$ for the minimisation was 19, having 6 degrees of freedom (the probability of obtaining a value for $`\chi ^2`$ bigger than 19 is $`0.0042`$). We further try to fit the data with two perturbers in the system. If we put two equal perturbers into the system (i.e. both with fixed $`\theta _\mathrm{c}=1\mathrm{mas}`$ and $`\sigma _{\mathrm{NIS}}=10\mathrm{km}/\mathrm{s}`$), we get a $`\chi ^2`$ of 15 (see Table 5) with 4 degrees of freedom. The probability of obtaining a $`\chi ^2`$ larger than that value is now $`0.0047`$. Such a reduction of $`\chi ^2`$ is apparently not a significant improvement of the fit (compared to the model with a single perturber).

Just for comparison, we also include the deconvolved image shapes in the fit with SIE+SH model. The resulting $`\chi ^2`$ is 140 for 8 degrees of freedom (the probability of obtaining a value for $`\chi ^2`$ larger than that is now only $`10^{28}`$). We essentially get the same model as when only fitting image positions and fluxes (see Table 6). We see that the fluxes provide much stronger constraints (their uncertainty is much smaller) than the image shapes.

Apparently, models with substructure yield significantly better fits than the ones without. However, we should stress that, strictly speaking, the appropriate statistical treatment can not be easily performed because the model parameters enter the $`\chi ^2`$ function in a non-linear way. As a result, the $`\chi ^2`$ function does not behave as a chi-square random variable with the appropriate number of degrees of freedom. One can see that already from the fact that the $`\chi ^2`$ function did not vanish when using a model with zero degrees of freedom. Hence, the probabilities we quoted above have to be taken with care; still they can be used to compare the performance of different models.

It is clear that it is difficult to simultaneously stretch and rotate the images with one (or two) perturber(s). The macro-model is not very successful in predicting both components of the image ellipticities and the fluxes, and therefore corrections are needed in the case of all four images. We can safely assume that the inclusion of further subclumps in the model would eventually lead to a “perfect” fit with the observed data. In particular three more subclumps close to the images B, C, and D would yield a significant improvement to the $`\chi ^2`$.

## 6 Strong lensing by an N-body simulated galaxy

A question that arises from model fitting of B1422+231 is whether such behaviour is seen by the N-body simulated galaxy and therefore generic of a typical galaxy lens. We used the cosmological N-body simulation data including gasdynamics and star formation of Steinmetz & Navarro (2001) for this purpose.

The simulations were performed using GRAPESPH, a code that combines the hardware N-body integrator GRAPE with the Smooth Particle Hydrodynamics technique (Steinmetz 1996). GRAPESPH is fully Lagrangian and optimally suited to study the formation of highly non-linear systems in a cosmological context. The version used here includes the self-gravity of gas, stars, and dark matter components, a three-dimensional treatment of the hydrodynamics of the gas, Compton and radiative cooling (assuming primordial abundances), the effects of a photo-ionizing UV background, and a simple recipe for transforming gas into stars.

The original simulated field is located at $`z=0.2`$ and contains approximately 300000 particles. The simulation is contained within a sphere of diameter $`32\mathrm{Mpc}`$ which is split into a high resolution sphere of diameter $`2.5\mathrm{Mpc}`$ centred around a galaxy and an outer low resolution shell. Gas dynamics and star formation is restricted to the high resolution sphere (280000 particles, 92000 of which dark matter), while the 34000 dark matter particles of the low resolution sphere sample the large scale matter distribution in order to appropriately reproduce the large scale tidal fields (see Navarro & Steinmetz 1997 and Steinmetz & Navarro 2000 for details on this simulation technique).

The simulation was performed in a $`\mathrm{\Lambda }`$CDM cosmology ($`\mathrm{\Omega }_0=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`\mathrm{\Omega }_\mathrm{b}=0.019/h^2`$, $`\sigma _8=0.9`$). It has a mass resolution of $`1.26\times 10^7M_{}`$ and a spatial resolution of $`0.5\mathrm{kpc}`$. A realistic resolution scale for an identified substructure is typically assumed to be $`40`$ particles which corresponds to $`5\times 10^8M_{}`$. The quoted mass resolution holds for gas/stars. The high resolution dark matter particles are about a factor of 6 ($`=\mathrm{\Omega }_0/\mathrm{\Omega }_\mathrm{b}`$) more massive.

From the original simulated field we took a cut-out map of size $`160\times \mathrm{\hspace{0.33em}160}\mathrm{kpc}`$ that is centred on a single galaxy. This area lies well within the high resolution sphere and is void of any massive intruder particles from the low resolution shell. The resulting mass distribution was smoothed using convolution with a Gaussian kernel characterised by a standard deviation of $`\sigma 0.8\mathrm{kpc}0.2\mathrm{arcsec}`$.<sup>1</sup><sup>1</sup>1For the distance calculations through the paper we assumed an Einstein-de-Sitter Universe and the Hubble constant $`H_0=65\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$. This map was then evaluated on $`2048\times \mathrm{\hspace{0.33em}2048}`$ grid points. The final map contains information about approximately 130000 particles with a total mass of $`3.0\times 10^{12}\mathrm{M}_{}`$. The surface mass density $`\kappa `$ was calculated for every grid point. We chose the redshift of the source to be $`z=3`$. A part of the cut-out map can be seen in Fig. 1.

The lens properties are calculated on a grid of $`2048\times \mathrm{\hspace{0.25em}2048}`$ points. The Poisson equation

$$^2\psi (𝜽)=2\kappa (𝜽).$$
(10)
is solved (on the grid) in Fourier space by the FFT (Fast Fourier Transformation) method. This method is incorporated in the routine KAPPA2STUFF from the IMCAT software of Nick Kaiser (http://www.ifa.hawaii.edu/~kaiser) that we used. It takes a grid map of $`\kappa (𝜽_j)`$ as an input and returns the values of deflection angle and complex shear (in x-space), along with other quantities. From these data we calculated the Jacobi matrix for each grid point.

The simulated galaxy is a field galaxy. Therefore we add two external shear components to the Jacobi matrix (evaluated at each grid point) in order to make it similar to the galaxy in B1422+231. The shear components were taken to be the same as the ones obtained from the best fitting SIE+SH model

$$\gamma _1^{\mathrm{ext}}=0.04,\gamma _2^{\mathrm{ext}}=0.16.$$

The external shear accounts for the effect of the neighbouring galaxies of the compact group, which are not present in the simulation.

Fig. 2 shows the magnification map of the surface mass density (given in Fig. 1) with additional external shear. One can clearly see the outer critical curve (white curve), while only the traces of the inner critical curve are visible (little circle inside the black region). The reason why we can see the outer critical curve so much better than the inner one is the following. At the centre $`𝜽_\mathrm{c}`$ of the galaxy, the surface mass density is very high and the determinant of the Jacobi matrix can be approximated by $`detA\kappa (𝜽_\mathrm{c})^2`$. In our particular case, $`\kappa (𝜽_\mathrm{c})50`$, and since at the critical curve $`detA=0`$, we see that the determinant has to decrease from 2500 to 0 in a region of $`0.4\mathrm{arcsec}`$. The transition is therefore very steep and we have a very good chance to miss the maximum value of magnification there, since the resolution is not high enough. At the outer critical curve, the change is slower and we can clearly see the points of high magnification.

In order to generate a similar image configuration as the one in B1422+231, one considers the caustic curve. This can be done by simply mapping the points of high magnification onto the source plane. Such a map is presented in Fig. 3.

On first sight the caustic structure of the N-body simulated galaxy looks the same as e.g. the caustic of the smooth NIE model with a small core radius. However, if we look only at the inner, asteroid caustic we can see that it is not completely “smooth”. With the help of bilinear interpolation we recalculated the magnification map on a refined grid (increasing the number of points being evaluated by $`5\times 5`$) and the corresponding caustic for such a grid is shown in Fig. 4. The caustic structure is much more complicated than in the case of a smooth model; in addition to folds and cusps we also have swallowtails formed (see Schneider et al. 1992 for more details on catastrophe theory).

### 6.1 Generating synthetic data

We select a source position $`𝜷_\mathrm{s}`$ such as to lie inside the asteroid caustic and close to the cusp, trying to chose a position for which we would get similar flux ratios as in the case of B1422+231. In total we considered 11 different source positions (see Fig. 4).

For each of them we first determine approximate image positions using the method described in Sect. 5.1. In order to get exact image positions we use the root finding method MNEWT again, for which we need the deflection angle to be continuous inside the region where we look for images. In our case the deflection angle is defined only on a $`2048\times 2048`$ grid. We perform bilinear interpolation between grid points. Having the image positions, we performed bilinear interpolation of the magnification map in order to get more exact magnification factors.

### 6.2 Fitting the synthetic data

For the $`\chi ^2`$-fitting method according to (7) we need to determine the uncertainty on the image positions and fluxes. Since we use interpolation for the MNEWT method we do not have a real estimate for the errors. One can, for example, set the errors to the same (relative) values as the uncertainties on the observed radio positions in B1422+231. For the typical scales we are using here (i.e. the distance B to D is approximately $`3.5\mathrm{arcsec}`$) this would mean an uncertainty of much less than a distance between two grid points. However, such a small error estimate is not realistic; due to the finite grid we estimate the image position uncertainties to be the distance between two grid points (which is a generous estimate; we use bilinear interpolation so the uncertainty is probably lower). The flux ratio errors were then set to be approximately 2 % for images A, B and C and 5 % for image D. These uncertainties are set to be the same as in the case of observed radio fluxes in B1422+231. The galaxy position error is set twice as big as that for the image positions.

The fitting procedure is performed in the same way as for the B1422+231 data. Again, image B is taken as a reference and the $`\chi ^2`$-function is evaluated according to (7). We try to fit the positions and fluxes with SIE+SH and SIE+SH+NIS models. The flux ratios of the 11 sets of synthetic data, together with the results of the minimisation are presented in table 7. We experience similar problems fitting fluxes as before; the $`\chi ^2`$-function value is high for all 11 data sets.

What might be surprising is the fact that we do not recover some properties of the lensing system. We know that the primary lens is fairly circular (one can not see that from fig. 2, since there, external shear is already added) and we know the values of the external shear components. What we also know a priori are the magnification factors for the images. These values were not recovered with high accuracy in model fitting.

For a smooth model and a source close to the cusp, the flux relation described before holds. We see that the fluxes violate this rule in all configurations we used. As we mentioned before this relation can only be violated when the source is away from the cusp or if there is substructure in the system. Since here we know the source position, the N-body lensing results show that the substructure we have in this particular simulated galaxy is indeed responsible for the observed deviation.

### 6.3 Discussion of the results from N-body lensing

An important question is whether the N-body simulation galaxy we are using is a good representation of a real galaxy for the purpose of lensing. If the resolution is not high enough, an N-body simulated galaxy might show more substructure than a real galaxy has.

In order to obtain the surface mass density map representative of lensing and to try to make sure that the substructure we see is not of numerical origin we used a smoothing length for particles of $`\sigma =0.8\mathrm{kpc}0.2\mathrm{arcsec}`$. The individual mass clumps we see in the corresponding surface mass density map (Fig. 1) therefore contain well above 100 particles. As we mentioned before, a realistic resolution scale for an identified substructure in the simulation corresponds to $`40`$ particles.

The regions that are interesting for multiple image formation typically have $`\kappa `$ values of order unity. In such regions we have approximately 300 particles per smoothing area, which, if the noise were Poissonian, would result in an error of about $`1/\sqrt{N}0.06`$. However, studies of errors of N-body simulations using smoothed particle hydrodynamics (Monaghan 1992; Niedereiter 1978) show that the errors are significantly smaller, and behave as $`\mathrm{log}N/N`$.

Any deviation of $`\kappa `$ due to substructure larger than the estimated error of the N-body simulations tells us that we are dealing with “physical” substructure, i.e. not of numerical origin. Our surface mass density maps show deviations well above the estimated Poisson noise, and therefore also above numerical noise. The results shown in Table 7 are affected by the “physical” substructure as well as by the noise. Although the expected noise is significantly smaller than believable substructure observed in the N-body simulations, it is hard to precisely estimate the relative contributions in the results shown here.

There are indications that the level of substructure as obtained from simulations can influence lensing phenomena a lot. In particular, the synthetic fluxes we obtained deviate strongly from those predicted by smooth models. This particular example of a simulated galaxy can of course not give us the answers to the aforementioned questions. In order to draw stronger conclusions, one would have to investigate many different realisations of N-body simulated galaxies and in addition use higher resolution simulations (currently unavailable). A statistical analysis to investigate the strong lensing properties could then be made. This is, however, beyond the scope of this paper and will be considered in a future work.

## 7 Conclusion

In this work we have investigated the influence of substructure in the gravitationally lensed system B1422+231. While it is intuitively clear that a lens galaxy is not a smooth entity, we have tried to investigate how deviation from a smooth model can influence lensing phenomena, especially the image flux ratios.

We have used two different smooth models for the lensing galaxy (SIE+SH and NIE+SH), and both failed very badly in fitting the image fluxes (we got $`\chi ^2=130`$ with 2 degrees of freedom). The use of models with substructure requires additional observational constraints. Therefore, we used deconvolved image shapes as constraints. We get a significant improvement of the fit compared to the smooth model. However, the way the substructure is introduced is oversimplified, thus we should not be surprised that the resulting $`\chi ^2`$ is still high. For the model with a single perturber we got $`\chi ^2=19`$ for 6 degrees of freedom, and with two perturbers we had $`\chi ^2=15`$ for 4, while the model without substructure (where deconvolved image shapes were included) gives $`\chi ^2=140`$ for 8 degrees of freedom.

Up to now we have not considered the possibility that microlensing plays a role for the radio fluxes. Koopmans & de Bruyn (2000) claim that they have detected microlensing in the multiply-imaged radio source B1600+434. Microlensing is a very tempting explanation for difficulties in fitting the fluxes for it can also explain why the $`8.4\mathrm{GHz}`$ A:B flux ratio has changed from 0.97 in 1991 (Patnaik et al. 1992) to 0.93 in 1997 (Patnaik et al. 1999). This has a consequence that again speaks in favour of substructure, since the presence of radio microlensing indicates that there is a significant number of compact objects in the lens galaxy halo.

N-body simulation data of a model galaxy provides a test for the influence of mass-substructure in strong gravitational lensing. When we generated data of four image systems with the simulated galaxy we again experienced difficulties. We have tried to fit image positions and fluxes and failed to obtaining a model that fits well. From these experiments we can conclude that the level of substructure obtained from this particular N-body simulated galaxy can cause the same difficulties as experienced in some of the real gravitationally lensed systems.

In order to obtain stronger conclusions one would have to investigate more realisations of simulated galaxies, also at different redshifts. However, the N-body simulation results indicate that substructure plays an important role in strong lensing. Also, modelling B1422+231 shows that the fluxes of more than one image are probably affected by the lens substructure. With the conclusions one can draw from lens modelling we can not give a direct answer about the properties of mass substructure in B1422+231.

One should therefore avoid using the flux constraints directly; they should, rather, be treated in statistical manner, e.g. in a way suggested by Mao & Schneider (1998). Fortunately, the perturbations on the scales we are dealing with here do not influence the image positions significantly, and play even less of a role for the time delay (Mao & Schneider 1998). Strong lensing thus remains one of the best tools to constrain the Hubble constant.

###### Acknowledgements.

We would like to thank the referee for many constructive comments which have helped to greatly improve this paper. Further we would like to thank Douglas Clowe and Alok Patnaik for many useful discussions. This work was supported by the Bonn International Physics Programme, by the Deutsche Forschungsgemeinschaft, and by the TMR Network “Gravitational Lensing: New Constraints on Cosmology and the Distribution of Dark Matter” of the EC under contract No. ERBFMRX-CT97-0172.