text:
# Linearized perturbation on stationary inflow solutions in an inviscid and thin accretion disc

## 1 Introduction

The object of this work has been to study the stability of the stationary inflow solutions in axially symmetric accretion, under the influence of a linearized perturbation. The chosen model has been that of an inviscid and thin accretion disc.

Use of this model, in which viscous effects have been ignored, has found a fairly widespread and consistent favour over the years (Abramowicz & Zurek, 1981; Chakrabarti, 1989; Molteni et al., 1996). As a consequence of such a physical prescription, there remains no mechanism for the outward transport of the angular momentum of the flow, and hence the specific angular momentum of the flow has been taken to be a constant of the motion. This naturally simplifies the analysis to a great extent, without significantly detracting from a true understanding of the underlying physics. The whole analysis has been based on the thin-disc configuration in accretion studies – a very well-known textbook model (Pringle, 1981; Frank et al., 1992) that has been quite readily and regularly invoked by researchers in accretion astrophysics.

It has been argued here, that in an inviscid and thin accretion disc, the stationary inflow solutions of abiding interest, are all stable under the influence of a linearized time-dependent perturbation. The technique adopted has been similar to the one already existing for the case of spherically symmetric accretion (Petterson et al., 1980; Theuns & David, 1992). Separate emphasis has been laid on discussing the stability of the so-called transonic inflow solution – a solution that passes through the critical point(s) of the flow, and whose velocity is subsonic very far away from the accretor, becoming supersonic as the accretor is approached. The subsonic inflows, a class of physically meaningful inflow solutions, the maximum of whose velocity always remains below the velocity of the inflow solution passing through the critical point(s), have been quite extensively considered as regards their stability, since, like the transonic inflow, the subsonic inflow solutions also obey the same outer boundary condition that the flow velocity should decrease very far away from the accretor. And it may also be pointed out that since the subsonic solutions do not pass through any singular point of the flow, the mathematical problem of studying their stability is much simpler compared with that of the transonic solution.

The classical Newtonian potential has been taken as the driving factor which effects the infall of matter, although it is easy to see that the use of the pseudo-Newtonian potential of Paczyński & Wiita (1980) – which satisfactorily represents general relativistic effects in a Newtonian framework and which is frequently in use to study the infall of matter under the gravitational influence of a black hole – alters nothing radically about the conclusions drawn from a linear stability analysis of the inflow solutions. What it changes, however, is the number of critical points of the flow. In the context of the use of the Paczyński-Wiita potential, mention may be made here of the fact that the number of critical points of the flow is also determined by various constants of the motion like accretion rate, specific angular momentum and the specific energy of the flow. It has been shown (Chakrabarti, 1990) that for a given angular momentum of the flow, there exists a range of accretion rate, or equivalently a range of energy, such that the flow has three critical points. Of these, the two outer points are saddle points, flanking between themselves, a third centre-type point. The flow, however, has only one critical point for values of accretion rate and energy which do not lie within the given range that will develop three critical points.

In solutions which admit of more than one saddle-type critical point, the existence of shocks is a possibility (Chakrabarti, 1990). For accretion on to a black hole, a solution would then pass through the outermost critical point (which is a saddle point), undergo a shock and then pass through the innermost critical point (again a saddle point) to cross the horizon with supersonic velocities. For accretion on to a neutron star, a flow passing through the inner critical point must undergo another shock to meet the inner boundary condition. The behaviour of these shocks are understood by invoking the Rankine-Hugoniot conditions (Chakrabarti, 1989). Once again it is to be emphasized that the outer boundary condition for all kinds of physically meaningful inflow solutions – shocked or continuous – is that the bulk velocity should decrease to zero very far away from the accretor.

## 2 The Equations of the Flow

For the radial drift velocity $`v`$ and density $`\rho `$, the flow is governed by Euler’s equation and the continuity equation. The radial component of Euler’s equation is given by

$$\frac{v}{t}+v\frac{v}{R}+\frac{1}{\rho }\frac{P}{R}+\frac{V}{R}\frac{L^2}{R^3}=0$$
(1)
where $`L`$, which is a constant, is the specific angular momentum, and the Newtonian gravity potential $`V=GM/R`$. The choice of the pseudo-Newtonian potential of Paczyński & Wiita (1980) would give $`V=GM/(Rr_g)`$, in which $`r_g`$ is the Schwarzschild radius of a black hole.

Use has been made of a polytropic equation of state

$$P=k\rho ^\gamma $$
(2)
where $`\gamma `$ is the polytropic exponent, whose admissible range $`(1<\gamma <5/3)`$ is restricted by the isothermal limit and the adiabatic limit respectively.

A surface density $`\mathrm{\Sigma }`$ has been defined by vertically integrating $`\rho `$ over the disc thickness $`H`$, and in the thin-disc approximation (Frank et al., 1992) this gives $`\mathrm{\Sigma }\rho H`$. The continuity equation, in terms of the surface density $`\mathrm{\Sigma }`$, in the thin-disc approximation, is thus given by

$$\frac{\mathrm{\Sigma }}{t}+\frac{1}{R}\frac{}{R}\left(\mathrm{\Sigma }vR\right)=0$$
(3)
Again invoking the thin-disc approximation, the disc thickness $`H`$ is obtained from the vertical component of Euler’s equation as

$$H\frac{c_s}{(GM)^{1/2}}R^{3/2}$$
(4)
where $`c_s`$, the speed of sound, obeys the relation

$$c_{s}^{}{}_{}{}^{2}=\gamma k\rho ^{\gamma 1}.$$
(5)
The set of equations (1)–(5), defines the whole problem completely. From (1) and (2) is obtained

$$\frac{v}{t}+v\frac{v}{R}+k\gamma \rho ^{\gamma 2}\frac{\rho }{R}+\frac{V}{R}\frac{L^2}{R^3}=0$$
(6)
while (3),(4) and (5) give

$$\frac{}{t}\left[\rho ^{(\gamma +1)/2}\right]+R^{5/2}\frac{}{R}\left[\rho ^{(\gamma +1)/2}vR^{5/2}\right]=0$$
(7)
The two equations (6) and (7) above are those that govern the flow in the inviscid and thin accretion disc.

## 3 The Perturbation Equation

The steady-state solutions of (6) and (7) are $`v_0`$ and $`\rho _0`$, on which are imposed small first-order perturbations $`v^{}`$ and $`\rho ^{}`$ respectively. Following a similar prescription by Petterson et al. (1980) for spherically symmetric accretion, it is convenient to introduce a new variable, $`f\rho ^{(\gamma +1)/2}vR^{5/2}`$, which in physical terms is the mass flux. Its steady-state value $`f_0`$, as can be seen from (7), is a constant that is identified as the mass accretion rate. A linearized perturbation $`f^{}`$, about $`f_0`$ is given by,

$$f^{}=\rho _{0}^{}{}_{}{}^{(\gamma 1)/2}R^{5/2}\left(v^{}\rho _0+\frac{\gamma +1}{2}v_0\rho ^{}\right)$$
(8)
Linearizing in the perturbation variables $`v^{}`$ and $`\rho ^{}`$, gives from (7), the expression,

$$\frac{\rho ^{}}{t}+\frac{2}{\gamma +1}\rho _{0}^{}{}_{}{}^{(\gamma 1)/2}R^{5/2}\frac{f^{}}{R}=0$$
(9)
With the use of (9), successive differentiation of (8) with respect to time yields, first

$$\frac{v^{}}{t}=R^{5/2}\rho _{0}^{}{}_{}{}^{(\gamma +1)/2}\left(\frac{f^{}}{t}+v_0\frac{f^{}}{R}\right)$$
(10)
and then

$$\frac{^2v^{}}{t^2}=R^{5/2}\rho _{0}^{}{}_{}{}^{(\gamma +1)/2}\left[\frac{^2f^{}}{t^2}+v_0\frac{}{R}\left(\frac{f^{}}{t}\right)\right]$$
(11)
Now linearizing in terms of the perturbation variables, gives from (6) the equation,

$$\frac{v^{}}{t}+\frac{}{R}\left(v_0v^{}+c_{s0}^{}{}_{}{}^{2}\frac{\rho ^{}}{\rho _0}\right)=0$$
(12)
where $`c_{s0}`$ is the speed of sound in the steady state.

Differentiating (12) with respect to time, and using the conditions given by (9),(10) and (11) will finally deliver the perturbation equation as

$$\frac{^2f^{}}{t^2}+2\frac{}{R}\left(v_0\frac{f^{}}{t}\right)+\frac{1}{v_0}\frac{}{R}\left[v_0\left(v_{0}^{}{}_{}{}^{2}\frac{2}{\gamma +1}c_{s0}^{}{}_{}{}^{2}\right)\frac{f^{}}{R}\right]=0$$
(13)
This expression bears a very close resemblance in form to a perturbation equation obtained by Theuns & David (1992), using similar methods, for the case of spherically symmetric accretion.

## 4 Stability Analysis of Solutions

Before taking up the stability analysis of solutions, as implied by (13), it would be worthwhile to have an understanding of the critical points of the flow. In the steady state, (6) and (7) can be combined to give the critical points – the required condition being that the numerator and the denominator of $`dv_0/dR`$ vanish simultaneously at those points. This leads to the expressions

$`v_{0}^{}{}_{}{}^{2}`$ $`=`$ $`{\displaystyle \frac{2}{\gamma +1}}c_{s0}^{}{}_{}{}^{2}`$
$`{\displaystyle \frac{5}{\gamma +1}}c_{s0}^{}{}_{}{}^{2}`$ $`=`$ $`{\displaystyle \frac{GM}{R}}{\displaystyle \frac{L^2}{R^2}}`$ (14)
The second equation above, being a quadratic in $`R`$, indicates that there would be two critical points for the flow. The choice of the pseudo-Newtonian potential of Paczyński & Wiita (1980), merely changes the number of critical points (Chakrabarti, 1989), without altering the first of the two equations above.

It is seen that the transonic inflow which originates very far away from the accretor, with a very low subsonic velocity, will pass through the first point – which is a saddle point, and then tend to curl around the second point – which is a centre-type point. In case the accretor is a large and distended object like an ordinary star (which justifies the use of the conventional Newtonian potential in this analysis), it is quite likely then that while the flow solution tends to curl around the inner critical point, it would also meet the surface of the accretor. This, as will be discussed later, has important implications for the stability of the transonic inflow.

Returning now to (13), a solution of the form $`f^{}=g(R)e^{\mathrm{i}\mathrm{\Omega }t}`$ is chosen, in which $`g(R)`$ is the spatial part of the perturbation and $`\mathrm{\Omega }`$ is real. This leads to the result

$$g\mathrm{\Omega }^22\mathrm{i}\mathrm{\Omega }\frac{d}{dR}(v_0g)+\frac{1}{v_0}\frac{d}{dR}\left[v_0\left(v_{0}^{}{}_{}{}^{2}\beta ^2c_{s0}^{}{}_{}{}^{2}\right)\frac{dg}{dR}\right]=0$$
(15)
where $`\beta ^2=2/(\gamma +1)`$.

For subsonic flows, it is easy to understand that there are two values of $`R`$ – one close to the accretor and one very far away – where the perturbation, which is in the form of a standing wave, can be constrained to die out. Multiplying (15) by $`v_0g`$ and integrating over the range of $`R`$, which is bounded by the two points where the perturbation dies out, yields,

$$\mathrm{\Omega }^2v_0g^2𝑑R=v_0\left(v_{0}^{}{}_{}{}^{2}\beta ^2c_{s0}^{}{}_{}{}^{2}\right)\left(\frac{dg}{dR}\right)^2𝑑R$$
(16)
The above result could be obtained because the boundary terms all vanish. With a real value of $`g`$, for the subsonic flows – all of them governed everywhere by the condition $`v_{0}^{}{}_{}{}^{2}<\beta ^2c_{s0}^{}{}_{}{}^{2}`$$`\mathrm{\Omega }^2`$ will be positive, implying that the perturbation, which is in the nature of a standing wave, will be oscillatory and at least will not grow in time. This line of reasoning was analogously established for spherically symmetric subsonic flows by Petterson et al. (1980).

It is also possible to subject the subsonic flow solutions to an analysis where the perturbation may be treated as a travelling wave. The method for this, again follows a prescription used by Petterson et al. (1980) for the spherically symmetric flow. The perturbation would have to be confined to a short-wavelength regime, in which the wavelength of the propagating waves will be much less compared to a characteristic length-scale of the system, which, in this case, as Petterson et al. (1980) argued for the analogous spherically symmetric situation, could be the radius of the accretor, and that for an ordinary star is $`R_{}`$. Since such a length-scale would be significantly large, it is somewhat justifiable to restrict the study to a perturbation of comparatively short wavelengths. Consequently, $`\mathrm{\Omega }`$, the frequency of the waves, would be large.

It is convenient to recast (15) as

$$\left(v_{0}^{}{}_{}{}^{2}\beta ^2c_{s0}^{}{}_{}{}^{2}\right)\frac{d^2g}{dR^2}+\left[3v_0\frac{dv_0}{dR}\beta ^2\frac{1}{v_0}\frac{d}{dR}\left(v_0c_{s0}^{}{}_{}{}^{2}\right)2\mathrm{i}\mathrm{\Omega }v_0\right]\frac{dg}{dR}\left(2\mathrm{i}\mathrm{\Omega }\frac{dv_0}{dR}+\mathrm{\Omega }^2\right)g=0$$
(17)
A solution to (17) may be obtained by expanding $`g(R)`$ as a power series (Petterson et al., 1980), in the form

$$g(R)=\mathrm{exp}\left[\underset{n=1}{\overset{\mathrm{}}{}}\frac{k_n(R)}{\mathrm{\Omega }^n}\right]$$
(18)
and to solve for the first few terms.

For such a solution, terms in $`\mathrm{\Omega }^2`$, $`\mathrm{\Omega }`$ and $`\mathrm{\Omega }^0`$ will respectively yield the conditions

$$\left(v_{0}^{}{}_{}{}^{2}\beta ^2c_{s0}^{}{}_{}{}^{2}\right)\left(\frac{dk_1}{dR}\right)^22\mathrm{i}v_0\frac{dk_1}{dR}1=0$$
(19)
$$\left(v_{0}^{}{}_{}{}^{2}\beta ^2c_{s0}^{}{}_{}{}^{2}\right)\left(\frac{d^2k_1}{dR^2}+2\frac{dk_1}{dR}\frac{dk_0}{dR}\right)+\left[3v_0\frac{dv_0}{dR}\frac{\beta ^2}{v_0}\frac{d}{dR}\left(v_0c_{s0}^{}{}_{}{}^{2}\right)\right]\frac{dk_1}{dR}2\mathrm{i}v_0\frac{dk_0}{dR}2\mathrm{i}\frac{dv_0}{dR}=0$$
(20)
and

$$\left(v_{0}^{}{}_{}{}^{2}\beta ^2c_{s0}^{}{}_{}{}^{2}\right)\left[\frac{d^2k_0}{dR^2}+2\frac{dk_1}{dR}\frac{dk_1}{dR}+\left(\frac{dk_0}{dR}\right)^2\right]+\left[3v_0\frac{dv_0}{dR}\frac{\beta ^2}{v_0}\frac{d}{dR}\left(v_0c_{s0}^{}{}_{}{}^{2}\right)\right]\frac{dk_0}{dR}2\mathrm{i}v_0\frac{dk_1}{dR}=0$$
(21)
It is to be noted that as compared to the spherically symmetric case, the analogous results for the inviscid and thin accretion disc, given by (19), (20) and (21), vary by a scale factor of $`\beta `$ for the speed of sound in the steady state. It is then quite straightforward to arrive at the expressions,

$$\frac{dk_1}{dR}=\frac{\mathrm{i}}{v_0\pm \beta c_{s0}}$$
(22)
and

$$k_0=\frac{1}{2}\mathrm{ln}\left(\beta v_0c_{s0}\right)+\mathrm{constant}$$
(23)
The first two terms in the power series expansion of $`g(R)`$ are given by (22) and (23). As with the spherically symmetric case, here also self-consistency may be had by establishing

$$\mathrm{\Omega }|k_1(R)||k_0|\mathrm{\Omega }^1|k_1(R)|$$
(24)
From the way it has been chosen, $`\mathrm{\Omega }(v_0\pm \beta c_{s0})/R_{}`$, it can be seen that the condition given by (24) is satisfied, if the radial distance is not too small. For large $`R`$, the asymptotic behaviour of $`k_1`$, $`k_1`$ and $`k_0`$ obtained from (21), (22) and (23) respectively, is given by

$$k_1\frac{1}{R},k_1R,k_0\mathrm{ln}R$$
(25)
The result given above upholds the self-consistency required by the condition in (24). The results in (25) follow from a very possible power-law dependence of the steady flow velocity on the radius, at great distances from the accretor.

The expression for the perturbation may then be written as

$$f^{}=\frac{e^{\mathrm{i}\mathrm{\Omega }t}}{\sqrt{\beta v_0c_{s0}}}\left[A_+\mathrm{exp}\left(\mathrm{i}\mathrm{\Omega }\frac{dR}{v_0+\beta c_{s0}}\right)+A_{}\mathrm{exp}\left(\mathrm{i}\mathrm{\Omega }\frac{dR}{v_0\beta c_{s0}}\right)\right]$$
(26)
In (26) there is a linear superposition of two solutions with arbitrary constants $`A_+`$ and $`A_{}`$, both travelling with velocity $`\beta c_{s0}`$, relative to the fluid, which itself moves with a bulk velocity $`v_0`$. One solution moves along with the bulk flow, while the other moves against it.

Use of (9) and (26) gives the perturbation in density as

$$\rho ^{}=\beta ^2R^{5/2}\rho _{0}^{}{}_{}{}^{(\gamma 1)/2}\frac{e^{\mathrm{i}\mathrm{\Omega }t}}{\sqrt{\beta v_0c_{s0}}}\left[\frac{A_+}{v_0+\beta c_{s0}}\mathrm{exp}\left(\mathrm{i}\mathrm{\Omega }\frac{dR}{v_0+\beta c_{s0}}\right)+\frac{A_{}}{v_0\beta c_{s0}}\mathrm{exp}\left(\mathrm{i}\mathrm{\Omega }\frac{dR}{v_0\beta c_{s0}}\right)\right]$$
(27)
Using (27), the perturbation in velocity can then be obtained from (8), and following some simple algebraic manipulation it can be rendered as

$$v^{}=\pm \frac{c_{s0}}{\beta }\frac{\rho ^{}}{\rho _0}$$
(28)
in which the positive and negative signs indicate outgoing and incoming waves, respectively.

In a unit volume of the fluid, the kinetic energy content is

$$_{\mathrm{kin}}=\frac{1}{2}\left(\rho _0+\rho ^{}\right)\left(v_0+v^{}\right)^2$$
(29)
The potential energy per unit volume of the fluid is the sum of the gravitational energy, the rotational energy and the internal energy, and is given by

$$_{\mathrm{pot}}=\left(\rho _0+\rho ^{}\right)\frac{GM}{R}\left(\rho _0+\rho ^{}\right)\frac{L^2}{2R^2}+\rho _0ϵ+\rho ^{}\frac{}{\rho _0}\left(\rho _0ϵ\right)+\frac{1}{2}\rho _{}^{}{}_{}{}^{2}\frac{^2}{\rho _{0}^{}{}_{}{}^{2}}\left(\rho _0ϵ\right)$$
(30)
where $`ϵ`$ is the internal energy per unit mass (Landau & Lifshitz, 1987).

The first-order terms vanish on time averaging. In that case, the contribution to the total energy in the perturbation comes from the second-order terms, which is given by

$$_{\mathrm{pert}}=\frac{1}{2}\rho _0v_{}^{}{}_{}{}^{2}+v_0\rho ^{}v^{}+\frac{1}{2}\rho _{}^{}{}_{}{}^{2}\frac{^2}{\rho _{0}^{}{}_{}{}^{2}}\left(\rho _0ϵ\right)$$
(31)
For an adiabatic perturbation, such as the one being prescribed here, the condition $`ds=0`$ can be imposed on the thermodynamic relation, $`dϵ=Tds+\left(P/\rho ^2\right)d\rho `$, which will then give

$$\frac{^2}{\rho _{0}^{}{}_{}{}^{2}}\left(\rho _0ϵ\right)|_s=\frac{c_{s0}^{}{}_{}{}^{2}}{\rho _0}$$
(32)
Combining the results from (28),(31) and (32) will give

$$_{\mathrm{pert}}=\frac{c_{s0}}{\beta ^2}\frac{\rho _{}^{}{}_{}{}^{2}}{\rho _0}\left[\frac{c_{s0}}{2}\left(\beta ^2+1\right)\pm \beta v_0\right]$$
(33)
The expression for $`\rho ^{}`$ is to be obtained from (27), and upon time averaging, the result will be

$$_{\mathrm{pert}}=\frac{1}{2}\frac{\beta A^2}{f_0}\frac{R^{5/2}\rho _{0}^{}{}_{}{}^{(\gamma 1)/2}}{\left(v_0\pm \beta c_{s0}\right)^2}\left[\frac{c_{s0}}{2}\left(\beta ^2+1\right)\pm \beta v_0\right]$$
(34)
where $`f_0\rho _{0}^{}{}_{}{}^{(\gamma +1)/2}v_0R^{5/2}`$, is a constant that in physical terms is the steady-state value of the accretion rate. The factor $`1/2`$ in (34 ) derives from the averaging of $`\rho _{}^{}{}_{}{}^{2}`$.

The total energy flux in the perturbation is obtained by multiplying $`_{\mathrm{pert}}`$ by the propagation velocity $`(v_0\pm \beta c_{s0})`$ and then by integrating over the area of the cylindrical face of the accretion disc, which is $`2\pi RH`$. Here $`H`$ is to be substituted from (4) and (5), and together with the fact that $`HR`$ in the thin-disc approximation (Frank et al., 1992), an expression for the energy flux will be delivered as

$$=\frac{\pi A^2}{f_0}\left(\frac{\gamma k}{GM}\right)^{1/2}\beta ^2\left[\pm 1+\frac{1\beta ^2}{2\beta (\pm \beta )}\right]$$
(35)
where $`=v_0/c_{s0}`$, is the Mach number.

It is to be seen that for purely subsonic flows, $``$ always remains less than $`\beta `$, which is a condition that is easily understood from the first of the two critical point conditions, given by (4). This would then imply that the total energy in the wave remains finite as the wave propagates. Indeed, for very great radial distances, when $`0`$, the flux is less compared to what it is near the outer critical point of the flow. The perturbation does not manifest itself as a runaway increase in the energy of the system. The system would then remain stable under the influence of such a perturbation.

For the case of the transonic flow, a different line of reasoning would have to be adopted, since the condition $`v_{0}^{}{}_{}{}^{2}<\beta ^2c_{s0}^{}{}_{}{}^{2}`$ would not hold everywhere, with the flow passing through a singular point. The physically meaningful flow in this case will be subsonic beyond the saddle-type outer critical point, pass through it to attain a supersonic velocity and then tend to curl around the inner centre-type critical point, but not actually flow through it. Theuns & David (1992) have established the stability of the transonic solution in spherical symmetry by arguing that the stability of the subsonic region of the flow, with the critical point as one boundary, will depend on conditions at the far outer end of the flow. In so far as the thin accretion disc is dependent only on $`R`$ (the disc being inviscid and hence, the related equations are independent of any angle variable), the analogous argument of Theuns & David (1992) is expected to hold good to ensure the stability of the subsonic region of the transonic solution.

For the supersonic region, which is finite, the reasoning of Garlick (1979) that a disturbance in this region will be carried away in a finite time, will also ensure the stability of the flow in this region.

There is another important argument to uphold the stability of the transonic flow in a thin accretion disc. In the case of the spherically symmetric transonic flow, the velocity near the accretor has an unbridled growth, such that any perturbation on it cannot be constrained to die out by a conceivable physical mechanism (Petterson et al., 1980). In the inviscid and thin accretion disc, the transonic solution has no such runaway growth in the inner region of the flow, tending to curl around the inner critical point as it does. And while doing so, it may reach the surface of the accretor, if the accretor is as distended as an ordinary star. This can ensure that a perturbation on the solution in this region will not behave in an unconstrained manner.

A note of caution has to be sounded here. The discussion presented above, pertains only to flows which are continuous everywhere. However, for spherically symmetric flows with discontinuities like shocks, Theuns & David (1992) have argued that the expectation would also be one in favour of the stability of those flows. Under the present assumptions governing an inviscid and thin accretion disc, the mathematical problem becomes very similar to the case of spherically symmetric accretion. In that case the argument of Theuns & David (1992) for spherical symmetry could very likely be extended to the axially symmetric system being considered here. In any case, regardless of the geometry of the problem, it is a matter of general understanding that a disturbance in the supersonic region in a shocked solution, is transmitted completely through the shock and that in the subsonic region a disturbance is reflected at the shock, with a diminished reflected amplitude for the disturbance as compared with the incident amplitude (Theuns & David, 1992). More to the point, by a local stability analysis, it has been demonstrated (Chakrabarti, 1989) that Rankine-Hugoniot shocks in an inviscid and thin accretion disc are stable under the influence of a short-wavelength perturbation in real time.

It can then be safely established that a physically well-behaved and meaningful inflow solution is indeed stable under the influence of a linearized perturbation. This conclusion holds good for both the transonic inflow solution, as well as the entire class of subsonic inflow solutions.

## 5 Concluding Remarks

So far, it has been demonstrated that linear stability analysis, which is essentially a perturbative technique, offers no clue about the exclusive choice of a particular solution – indeed, to the extent that a perturbative technique can be relied upon, every solution (transonic or subsonic) seems just as much realizable as any other. It is then to be conjectured that the criterion for the selection of a particular solution would have to be non-perturbative in character and would have to be based on fundamentally physical arguments, like the maximization of the accretion rate or the minimization of the total energy associated with a solution, as in the case of the transonic inflow solution in spherically symmetric accretion (Bondi, 1952; Garlick, 1979).

A recent work (Ray & Bhattacharjee, 2002) on spherically symmetric accretion has indicated that such an approach, however, would necessitate the addressing of the whole issue of the selection of a particular solution, from the viewpoint of the temporal evolution of the flow, instead of considering it solely on the basis of the stationary picture. The study has shown that it is the transonic solution which is decisively preferred to all the others. Such an expectation would also be valid for the case of the inviscid and thin accretion disc being discussed here.

Support for such a contention can also be had from a very different quarter. Study carried out over the last two decades has established that there is a very close one-to-one correspondence between certain features of black hole physics and the physics of supersonic acoustic flows. More specifically, for an inviscid, barotropic and irrotational fluid flow (such as the conventional spherically symmetric accretion) the equation describing an acoustic disturbance can be rendered in the form of a metric that relates very closely to the Schwarzschild metric, represented in the Painlevér-Gullstrand form (Visser, 1997). Close to the sonic point of the flow, the similarity is particularly evident. Such a close correspondence between the unidirectionality of trajectories near a black hole and of flows passing through the sonic point, is a reason strong enough to suggest that transonic trajectories are the favoured ones. Extending that argument, a comparison may also be drawn for rotational, inviscid fluid flows, to have an understanding of the favoured status of the trajectories in this case.

## Acknowledgments

This research has made use of NASA’s Astrophysics Data System. The financial assistance provided by the Council of Scientific and Industrial Research, Government of India, is being gratefully acknowledged. Gratitude is also to be expressed to Prof. J. K. Bhattacharjee for his very helpful advice, given unreservedly throughout this work.


format: text

text:
# Reinterpretation of the equilibrium magnetization of a Tl2Ba2CaCu2O8+x single crystal. Another phase transition in the mixed state of high-𝑇_𝑐 superconductors?

## Abstract

We apply a recently developed scaling procedure for the analysis of the equilibrium magnetization $`M`$ that was measured on a Tl2Ba2CaCu2O8+x single crystal and was recently reported in the literature. The results of our analysis are distinctly different from those obtained in the original publication where the Hao-Clem model served to analyze the magnetization data. We argue that the Hao-Clem model is not adequate for a quantitative description of the mixed state in high-$`T_c`$ superconductors especially in high magnetic fields. The scaled equilibrium magnetization data reveal a pronounced kink in the $`M(H)`$ dependence that might be indicative of a phase transition in the mixed state.

First, we briefly review the current situation in the interpretation of magnetization measurements in the mixed state of high-$`T_c`$ (HTSC) superconductors. The most widely used approach for evaluating the thermodynamic critical field $`H_c`$ and the Ginzburg-Landau (GL) parameter $`\kappa `$ of HTSC’s from magnetization data is the Hao-Clem model. 2 ; 3  Using the values of $`H_c`$ and $`\kappa `$ that result from the field dependence of the magnetization measured at different temperatures, the temperature dependence of the upper critical field $`H_{c2}(T)`$, as well as the magnetic field penetration depth $`\lambda (T)`$, and the superconducting coherence length $`\xi (T)`$ may be calculated. However a recent theoretical analysis of the Hao-Clem approach demonstrated that the calculation of the local magnetic induction by a linear superposition of contributions of individual vortices, as it was done in Ref. 3, , is quantitatively acceptable only for low magnetic fields. pog1 ; pog2  A rather direct confirmation of the inadequecy of the Hao-Clem model for the quantitative analysis of experimental results is provided by Figs. 3(a) and 3(b) in Ref. 4,  where data on the equilibrium magnetization of Tl-based single crystals were presented and analyzed. The quoted Fig. 3(a) reveals a rather strong increase of $`\kappa `$ with increasing temperature, resulting in a very unusual temperature dependence of the upper critical field with $`H_{c2}`$ being practically temperature independent at low temperatures and almost diverging at temperatures above 85-90 K (Fig. 3(b)). Because this anomalous increase of $`H_{c2}`$ starts well below the superconducting critical temperature $`T_c`$ of the investigated sample, it can hardly be ascribed to fluctuation effects. This very strong increase of $`\kappa `$ upon approaching the critical temperature is not a feature that was observed in one particular study. It turns out to be a rather universal result if the Hao-Clem model is employed for the analysis of magnetization data. hc1 ; hc2 ; hc3 ; hc4 ; hc5 ; hc6 ; hc7 ; hc8 ; hc9 ; hc10 ; hc11 ; hc12  We argue that the magnetic field dependence of the magnetization at constant temperature provided by the Hao-Clem model is not sufficiently accurate and therefore, a fitting of the relevant theoretical formula to experimental data results in an unphysical temperature dependence of the GL parameter $`\kappa `$, one of the key parameters in the Hao-Clem model. In the GL theory $`\kappa `$ is temperature independent and it only slowly decreases with increasing temperature in the framework of microscopic theories. th1 ; th2  Below we show that the quoted experimental data may very well be interpreted by assuming that $`\kappa `$ is temperature independent. The temperature dependence of $`H_{c2}`$ which results from our analysis, is in good agreement with expectations of the GL theory.

Our analysis of the original data which were presented in Ref. 4,  is based on a scaling procedure developed in Ref. 5, . It was shown that if the GL parameter $`\kappa `$ is temperature independent, the equilibrium magnetizations in the mixed state of type-II superconductors at two different temperatures $`T_0`$ and $`T`$ are related by

$$M(H/h_{c2},T_0)=M(H,T)/h_{c2}+c_0(T)H,$$
(1)
with $`h_{c2}(T)=H_{c2}(T)/H_{c2}(T_0)`$ representing a normalized upper critical field. 5  The first term on the right-hand side of Eq. (1) is universal for any type-II superconductor, while the second one is introduced to account for the temperature dependent paramagnetism of HTSC s in the normal state. 6  In the following we use $`M_{eff}`$ to denote the magnetization calculated using Eq. (1) in order to distinguish it from direct experimentally measured values of $`M`$. The quantities $`h_{c2}(T)`$ and $`c_0(T)`$ are the adjustable parameters of the scaling procedure and they are established from the condition that the $`M_{eff}`$ curves, calculated from the $`M(H)`$ measurements in the reversible regime at different temperatures, collapse onto a single curve, representing the equilibrium magnetization at $`T=T_0`$ (see Ref. 5,  for details).

It has been shown that the scaling procedure described by Eq. (1) works rather well for different HTSC materials including single crystals and grain-aligned samples, 5  as well as polycrystals and ceramics. 8  An important advantage of this approach is that no particular formula for $`M(H)`$ has to be assumed a priori and therefore, this scaling procedure may be used for any type-II superconductor, independent of the pairing type, the absolute value of $`\kappa `$, the anisotropy of superconducting parameters, or the demagnetization factor of the sample.

As has briefly been outlined above (see Ref. 5,  for details), Eq. (1) serves to evaluate the temperature variation of the normalized upper critical field $`h_{c2}(T)`$. It turns out that $`h_{c2}(T)`$ for this Tl2212 sample varies perfectly linearly with $`T`$ above 90 K. This linearity allows for a quite accurate evaluation of the critical temperature $`T_c`$ by extrapolating the $`h_{c2}(T)`$ curve to $`h_{c2}=0`$. The inset of Fig. 1 demonstrates that the value of $`T_c`$ resulting from this procedure is consistent with the temperature dependence of the low-field magnetization of this sample. The normalized upper critical field $`h_{c2}`$ versus $`T/T_c`$ is shown in the mainframe of Fig. 1.

As has been shown in our previous studies,5 ; 8  all investigated HTSC’s may so far be divided into two groups. For each group, the $`h_{c2}(T/T_c)`$ curves for the individual HTSC compounds turn out to be identical but the as yet identified two universal curves are distinctly different. The major group includes practically all different families of HTSC s, except Y2Ba4Cu7O15 and underdoped YBa2Cu3O7-x (YBCO) compounds which belong to the other group. As may be seen in Fig. 1, $`h_{c2}(T)`$ of this particular Tl2212 sample is practically identical with the $`h_{c2}(T)`$ curves of two underdoped YBCO samples investigated in Ref. 10, . This implies that the smaller group is not limited to underdoped Y-based cuprates. It was shown in Ref. 8,  that a ceramic Tl2212 sample with $`T_c=107`$ K, which was investigated in Ref. 11, , belongs to the other group represented by the dotted line in Fig. 1. A possible explanation for this difference between these two Tl-based samples may be sought in the difference of the oxygen content.

The dependence of $`M_{eff}(T_0)`$ on $`H/h_{c2}(T)`$ that results from our scaling procedure is shown in Fig. 2. Because of the high quality of the experimental data presented in Ref. 4,  and the extended covered range of magnetic fields, the scaling is nearly perfect. The $`M_{eff}(H)`$ data points, calculated from the measurements at different temperatures, combine to a single curve with virtually no deviations. The remarkable feature of this curve is a pronounced kink at $`H/h_{c2}20`$ kOe. This kink clearly indicates a significant change in the properties of the mixed state. Unfortunately, the measurements in Ref. 4,  are limited to magnetic fields $`H20`$ kOe and, as may be seen in Fig. 2, only a limited number of data points, measured at $`T96`$ K and $`H=20`$ kOe combine to the $`M_{eff}(H)`$ curve above the kink. This is why on the basis of the available data, no definite conclusions concerning the implications of this observation can be made and further studies are needed to clarify the situation. If this feature will be confirmed by a more detailed study, the kink would definitely reflect some kind of a phase transition.

We conclude that, contrary to the claims of Ref. 4, , the magnetization data perfectly well be explained by assuming a temperature independent Ginzburg-Landau parameter $`\kappa `$ and a conventional temperature dependence of the upper critical field$`H_{c2}`$, as shown in Fig. 1.


text:
# ANOMALOUS X-RAY PULSARS AND SOFT GAMMA-RAY REPEATERS — THE CONNECTION WITH SUPERNOVA REMNANTS

## INTRODUCTION

The anomalous X-ray pulsars (AXPs) and soft gamma-ray repeaters (SGRs) both represent exotic populations of neutron stars, distinct from more typical core-collapse products such as radio pulsars or X-ray binaries. We have learnt a great deal about SGRs and AXPs through their timing, photometric and spectroscopic properties — the consensus from these studies is that the SGRs, and most likely the AXPs also, are “magnetars”, neutron stars with surface magnetic fields $`10^{14}10^{15}`$ G (see reviews by Hurley 2001 and Mereghetti et al. 2003, and contributions in these proceedings by Woods and by Kaspi).

However, direct observations of AXPs and SGRs have as yet provided little data on their distances, ages, space velocities, lifetimes, birthsites and progenitors, all of which provide crucial information on the demographics and evolution of these interesting sources. In this review, I discuss how we can obtain indirect information on these issues by associating AXPs and SGRs with supernova remnants, star clusters and molecular clouds.

## ASSOCIATIONS WITH SUPERNOVA REMNANTS

Associations between neutron stars and supernova remnants (SNRs) can provide unique information about core-collapse supernovae and their aftermath. Over the last 35 years, there have been many controversial and colourful discussions about various radio pulsars and their associations (or lack thereof) with SNRs. Consequently, a considerable amount of collective wisdom has been accumulated as to how to best judge a potential association. Important criteria for considering an association include:

* Evidence for interaction. If a pulsar and a SNR can be shown to be interacting (e.g. if a pulsar is in the process of penetrating the shell of a SNR; Shull et al. 1989), this provides a very strong case for a physical association.
* Consistent distances/ages. Pulsar distances are usually estimated by comparing the dispersion of pulses with a model for the Galactic electron distribution, while a pulsar’s characteristic age is $`\tau _c=P/2\dot{P}`$ (where $`P`$ is the pulsar’s spin-period and $`\dot{P}`$ is the period derivative). SNR distances are most often determined from H i absorption measurements, while SNR ages can be estimated from their spectral properties or from an association with a historical supernova. For a genuine association, agreement between the pulsar’s and SNR’s distances are expected, and agreement in their ages is encouraging also.
* Inferred transverse velocity. If we assume that a neutron star was born at the geometric centre of a coincident or nearby SNR, then the offset of the pulsar from the SNR’s centre can be combined with an age estimate for the system to calculate an inferred transverse velocity for the pulsar. If the association is real, this inferred velocity should fall within the observed distribution of neutron star velocities (see Arzoumanian et al. 2002).
* Proper motion. In some cases, the proper motion of a pulsar can be measured. In a genuine association, the projected proper motion vector should be directed away from the SNR’s centre.
* Probability of random alignment. In the absence of other information, a basic indicator as to the likelihood of an association is the probability that the pulsar and the SNR lie near each other on the sky simply by chance. A pulsar sitting just outside a SNR in a complicated region of the inner Galaxy will usually correspond to a spurious association, but a pulsar located at a SNR’s centre in an otherwise empty region of the sky makes for a compelling case.

These criteria cannot all always be applied, and ultimately the judgement as to whether a potential association meets these requirements is a subjective one. If one lesson can be learnt from the claims and counterclaims of associations made over the last few decades, it is that an association is false until proven otherwise! With all this experience in hand, we can turn our attention to AXPs and SGRs, and their potential association with SNRs. However, we quickly find that most of the above criteria cannot be applied to these populations:

* To date, there is no evidence for any interaction of an AXP or SGR with an associated SNR, nor is it even clear whether such an interaction would produce any observable consequences.
* Independent distance estimates to AXPs and SGRs generally come only from measuring the foreground absorption in their X-ray spectra, and are correspondingly highly uncertain. The characteristic age for an AXP or SGR is $`\tau _c=P/2\dot{P}`$ as for radio pulsars. However, this formula is applicable only if a neutron star’s spin-down is entirely due to steady spin-down resulting from magnetic dipole radiation, while several of the AXPs and SGRs show complicated timing behaviour not consistent with this assumption (e.g. Kaspi et al. 2001a; Woods et al. 2002). Thus there is currently no way to make an independent and accurate distance or age estimate for an AXP/SGR.
* The transverse velocity for an AXP or SGR can be inferred from a potential SNR association, just as for radio pulsars. However, since there are no observational data on the underlying velocity distribution of AXPs and SGRs, there is nothing to which we can compare these inferred velocity estimates.
* Proper motion measurements of AXPs or SGRs would be extremely useful in assessing SNR associations, just as for radio pulsars. However, AXPs and SGRs are generally distant objects, and cannot be observed with VLBI techniques because of their lack of observed radio emission. Thus to date no AXP/SGR proper motions have been measured, nor are there any constraining upper limits.

Thus in considering SNR associations with AXPs and SGRs, all we are usually left with is the probability of random alignment. While this is still a useful criterion, extreme caution must be applied. As an example, in Figure 1 we show the radio emission from $`50`$ deg2 of the inner Galaxy, and have marked the positions of the 20 catalogued SNRs and 66 catalogued radio pulsars which fall within this region. If positional coincidence was our only guide, we would conclude that there are about a half dozen pulsar/SNR associations within this region. Application of the above criteria shows that most of these associations are erroneous, and only at most two or three are genuine.

## AXPS in SUPERNOVA REMNANTS

### Data On AXP/SNR Associations

There are five known AXPs, and a further two AXP candidates (see Mereghetti et al. 2003 for a review). Of these seven sources, three have convincing associations, as discussed by Gaensler et al. (2001) and listed in Table 1. In all three associations, the position of the AXP results in a very low probability of random alignment with the SNR, the age of the SNR is $`\stackrel{<}{_{}}`$10 kyr, and the offset of the AXP from the SNR centre implies a moderate transverse velocity for the neutron star, $`\stackrel{<}{_{}}500`$ km s-1.

The environments of the remaining four sources are shown in Figure 2. Of these sources, Gaensler et al. (2001) have shown that RX J170849–400910 is near a possible SNR but is in a very complicated region, that 1E 1048.1–5937 is near the Carina nebula but also has no associated SNR, and that 4U 0142+61 has no nearby sources. More recently, Lamb et al. (2002) have identified an AXP candidate, CXOU J010042.8–721132, in the Small Magellanic Cloud (SMC). As showm in the lower-right panel of Figure 2, CXOU J010042.8–721132 is near the H ii region and SNR complex N 66, but again has no specific SNR with which it can be associated.

One might draw the conclusion that the four AXPs without SNR associations are older, or are somehow different, from the three sources embedded in young SNRs. However, a note of caution is that only about half of all young radio pulsars have clear SNR associations. For example, the Crab pulsar, while powering a spectacular synchrotron nebula, has no surrounding SNR representing the ejecta from the associated supernova explosion (e.g. Frail et al. 1995). The lack of SNRs around some young pulsars is usually explained in terms of expansion into a cavity in the ISM or other low density region, and a similar explanation may well apply to the AXPs which lack SNRs.

### Implied Properties of AXPs

We have argued above that about half the AXPs (or AXP candidates) can be convincingly associated with SNRs. The associated AXPs all are located at their SNR’s centres, have inferred transverse velocities $`\stackrel{<}{_{}}500`$ km s-1, and have ages less than 10 kyr. In all four regards, these properties are very similar to those seen for the youngest radio pulsars and the so-called “central compact objects” (CCOs; see Pavlov et al. 2002). Figure 3 visually demonstrates how all three populations are found to be centrally located in young SNRs, implying that the AXPs are a population of young neutron stars. For radio pulsars, we see many associations in which the SNR is somewhat older and the pulsar is significantly offset from the SNR centre (e.g. Migliazzo et al. 2002). The fact that we do not see such associations for AXPs argues for an observable lifetime for AXPs of $`10`$ kyr. We can correspondingly infer a Galactic birth-rate for AXPs of $`>1/2000`$ yr (Gaensler et al. 1999).

If the braking torque on AXPs is due to electromagnetic dipole radiation, then a very high magnetic field, $`\stackrel{>}{_{}}10^{14}`$ G, is needed to slow AXPs from their presumed birth spin periods ($`P<1`$ s) to their observed long periods ($`P10`$ s) over their short lifetimes. Associations with young SNRs are thus consistent with the magnetar model for AXPs. In order to account for the limited AXP lifetime inferred from these associations, one must invoke a process such as magnetic field decay or rapid surface cooling (Colpi et al. 2000; Duncan 2002), which would render AXPs unobservable at later times.

Some authors have proposed that AXPs are accreting systems (e.g. Mereghetti & Stella 1995). However, such models have problems accounting for the optical/infrared properties of AXPs (Hulleman et al. 2000; Kern & Martin 2002), as well as for their recently discovered X-ray bursts (Gavriil et al. 2002; Kaspi & Gavriil 2002). Unusual accreting scenarios are also required to account for the SNR associations, because of the high braking torque needed to slow the AXPs to long periods over their short lifetimes (e.g. Chatterjee et al. 2000).

For a given age, the size of a SNR in the Sedov phase depends on the kinetic energy of the initial explosion and the density of the ambient medium. Marsden et al. (2001) have taken the SNRs associated with AXPs, estimated a distance to and age for each, adopted a fixed initial explosion energy of $`10^{51}`$ erg, and have consequently determined the ambient density in each case. These authors find that the densities into which the SNRs containing AXPs are expanding are consistently $`\stackrel{>}{_{}}0.1`$ cm-3, much higher than the density $`0.003`$ cm-3 typically seen around radio pulsars. Marsden et al. (2001) conclude that a high ambient medium around a supernova explosion is responsible for forming an AXP rather than producing a more normal radio pulsar.

However, there are some important deficiencies in such an argument (see Gaensler et al. 2001 and Duncan 2002 for more detailed discussions). First, the ages and explosion energies of SNRs are typically uncertain by up to a factor of two, while their distances are uncertain to $`20`$%. The resulting estimate of the ambient density is thus uncertain by two orders of magnitude. More fundamentally, there is a serious selection effect inherent in such studies, in that AXPs have mostly been discovered in targeted observations of SNRs, while most radio pulsars are found in all-sky surveys. Since observable SNRs are only produced in dense regions of the interstellar medium (e.g. Kafatos et al. 1980), it is not surprising that the known AXPs are found in such regions also. If one only considers radio pulsars associated with SNRs, one finds ambient densities $`0.2`$ cm-3 (Frail et al. 1994), no different from that found for AXPs in SNRs by Marsden et al. (2001). Conversely, if one could carry out an all-sky survey for AXPs, one would most likely find that most were in low-density regions. From the available evidence, it thus seems that the environments of the supernovae which produce AXPs are no different from those of the supernovae which produce radio pulsars.

## SGRS IN SUPERNOVA REMNANTS

### Data On SGR/SNR Associations

There are four known SGRs, plus one SGR candidate (see Hurley 2001 for a review). For all sources, the case for an association with a SNR is considerably less convincing than for the AXP/SNR associations discussed above.

SGR 0526–66 is on the rim of SNR N49 in the Large Magellanic Cloud (LMC) (Figure 4). This SNR is about 5000 yr old, an age comparable to that seen for the SNRs associated with AXPs. However, an important difference between this potential association and those discussed above is that the neutron star is significantly offset from the SNR centre, implying a transverse space velocity $`1000`$ km s-1. Over the whole LMC, the probability of a random alignment between these sources is $`0.002`$, but the system lies in a complicated region on the edge of the LMC-4 superbubble, in which several other adjacent SNRs have been identified (Figure 4; Gaensler et al. 2001). Thus while it seems likely that SGR 0526–66 is in a region of recent supernova activity, it is difficult to build a case for a specific association with N49.

SGR 1806–20 was originally found to be at the centre of the catalogued SNR G10.0–0.3. However, it has since been realised that the SGR is offset by $`14^{\prime \prime }`$ from the core of G10.0–0.3 (Kaplan et al. 2000a), while the latter’s amorphous morphology and time-variable behaviour suggests that the classification of G10.0–0.3 is erroneous, it rather being a nebula powered by the luminous blue variable star at its centre (Gaensler et al. 2001). Infrared and millimetre observations have subsequently shown that SGR 1806–20 is associated with a massive star cluster and molecular cloud complex, all at a distance of 15 kpc (e.g. Fuchs et al. 1999).

SGR 1900+14 is located $`5^{}`$ outside the rim of SNR G42.8+0.6 (see Figure 5). This is a complicated region of the inner Galaxy, and the probability of a random alignment between these two sources is $`4`$% (Gaensler et al. 2001; Kaplan et al. 2002b). While this alone makes the SGR/SNR association far from compelling, a young ($`\tau _c=38`$ kyr) radio pulsar, PSR J1907+0918 immediately adjacent to the SGR could just as likely be associated with the SNR (Lorimer and Xilouris 2000), further increasing the probability that at least one of these two neutron stars is not associated with the SNR. As is the case for SGR 1806–20, SGR 1900+14 is possibly associated with a massive star cluster, at a distance $`>10`$ kpc (Vrba et al. 2000).

SGR 1627–41 is near the SNR G337.0–0.1. The SNR is probably less than 5000 yr old (Sarma et al. 1997), and the SGR’s offset from the SNR correspondingly implies a projected space velocity $`>`$1000 km s-1 for the SGR. Again, the region is a very complicated one, Figure 6 demonstrating that both the SGR and SNR are embedded in the extended CTB 33 complex, which also contains OH maser emission, molecular clouds and several H ii regions all at a distance of $`11`$ kpc (e.g. Corbel et al. 1999). The probability of a random alignment in this case is high, although it is reasonable to suppose that the SGR and SNR are related to the same general episode of star formation and consequent supernova activity.

SGR 1801–23 is a candidate SGR which has not yet been well-localised, its error box passing near or through several SNRs in a particularly dense region of the Galactic plane (Cline et al. 2000). Until this source is confirmed and its position better refined, it is premature to propose an association of this source with any other object.

### Implied Properties of SGRs

While three out of four of the confirmed SGRs are in the vicinity of SNRs, in all cases the SGR is on the rim of or outside the SNR, implying a transverse velocity $`\stackrel{>}{_{}}1000`$ km s-1. In the absence of any other supporting evidence, none of these associations are very compelling.

However, the fact that SGRs are always near SNRs, as well as possibly being associated with massive star clusters and molecular cloud complexes, makes it clear that SGRs are broadly associated with recent star formation and supernovae. This argues that SGRs are neutron stars with ages $`<100`$ kyr, consistent with their interpretation as magnetars.

Further inferences about the nature of SGRs depend on how one interprets their potential association with SNRs. If one takes the three claimed SGR/SNR associations to be genuine, this implies that SGRs are a high velocity population, distinct from AXPs or radio pulsars. It has been proposed that newborn magnetars might experience asymmetric neutrino recoil (Duncan and Thompson 1992), in which case these anomalously high space velocities are not unexpected. However, this would imply that AXPs and SGRs are separate populations, a conclusion at odds with the many other similarities seen between these two groups of sources (e.g. Kulkarni et al. 2003).

Alternatively, if one concludes that the SGR/SNR associations are spurious, one can interpret the SGRs as a population of neutron stars with ages $`50`$ kyr, whose SNRs have dissipated. The SGRs thus must correspond to an older population than do the AXPs. If the AXPs and SGRs represent evolutionary phases in the life of a magnetar, then this suggests an evolutionary sequence in which AXPs evolve into SGRs (Gaensler et al. 2001). A problem with such a proposal is that it predicts that the spin periods of SGRs should generally be longer than those seen for the AXPs, which is not observed. One possibility is that the AXPs undergo extended periods of spin-up, but such episodes are yet to be seen in long-term monitoring observations (e.g. Gavriil & Kaspi 2002).

A final possibility is that the SGRs and AXPs represent approximately coeval populations of magnetar, but with different progenitor properties. For example, if SGRs result from very massive progenitors, then it is expected that the resulting supernova explosions will always occur in extended low-density wind bubbles. A SNR will be unobservable while expanding into such a bubble, and then would fade rapidly after it collides with swept-up material at the bubble’s edge (Braun et al. 1989). Such a possibility would account for the lack of SNR associations for the SGRs, and would also explain why some of these sources are associated with molecular clouds and massive star clusters.

Recent evidence suggests that AXPs and SGRs are most likely not distinct populations: two AXPs have shown SGR-like bursts (Gavriil et al. 2002; Kaspi & Gavriil 2002), and one SGR has a quiescent spectrum which is very similar to that of an AXP (Kulkarni et al. 2003). However, it is perhaps premature to claim that AXPs and SGRs represent the same underlying group of objects, distinguished only by the manner in which they were discovered. It seems more reasonable to suppose that there is a continuous distribution of properties across the magnetar population: one can speculate that the most massive progenitors result in SGRs, the most highly magnetised neutron stars which are consequently most prone to burst and which are least likely to be associated with SNRs. Lower mass progenitors produce AXPs, which are lower field neutron stars, which burst less often but which are generally associated with SNRs.

## CONCLUSIONS

Until a few years ago, there was a great deal of confusion both as to the environmental properties of SGRs and AXPs, and as to what these data were telling us about the nature of these exotic compact objects. However, in the last few years a series of multiwavelength efforts have made the picture more clear.

From these data, we can conclude that some AXPs have convincing SNR associations, which allows us to conclude that AXPs are young ($`<10`$ kyr) neutron stars with relatively low space velocities. On the other hand, the SGRs which have been identified are probably not associated with SNRs, implying that they are either older ($`\stackrel{>}{_{}}50`$ kyr) neutron stars whose SNRs have faded, or that they result from massive progenitors in low-density bubbles. The former case is consistent with a scenario in which the AXPs evolve into SGRs, while in the latter case the AXPs and SGRs can be coeval, but might be distinguished by the mass of their progenitor stars.

We are still dealing with a small number of sources, and in several cases the picture is still not completely clear. A number of efforts are being undertaken to clarify these issues. Kaplan et al. are undertaking an extensive multiwavelength project to find new neutron stars in SNRs, Hurley et al. (2001) are measuring the proper motion of SGRs with Chandra to determine whether their implied high space velocities are genuine, and Gaensler et al. are using H i data to search for low density cavities around SGRs. These efforts should provide new clues as to the nature of these intriguing objects.

## ACKNOWLEDGEMENTS

I thank David Kaplan, Patrick Slane and Miller Goss for providing me with images of regions around AXPs and SGRs. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The Australia Telescope is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO. This research has made use of the NASA/ IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

Email address of Bryan Gaensler: bgaensler@cfa.harvard.edu


text:
# Evolution of quantum systems with a scaling type of time-dependent Hamiltonians

## Abstract

We introduce a new class of quantum models with time-dependent Hamiltonians of a special scaling form. By using a couple of time-dependent unitary transformations, the time evolution of these models is expressed in terms of related systems with time-independent Hamiltonians. The mapping of dynamics can be performed in any dimension, for an arbitrary number of interacting particles and for any type of the scaling interaction potential. The exact solvability of a “dual” time-independent Hamiltonian automatically means the exact solvability of the original problem with model time-dependence.

PACS numbers: 03.65.-w, 03.65.Fd, 03.65.Ge, 02.30.Tb

Keywords: Quantum motion, time-dependent Hamiltonian, wave functions, unitary transformations, exactly solvable models.

Dynamics of quantum systems, governed by time-dependent Hamiltonians, has attracted much of attention for a long time. A relevant progress has been made mainly in the study of one-dimensional time-dependent harmonic oscillators1-10 which have many applications in various areas of physics (see e.g. Refs. 11-13). Lewis1 and Lewis and Riesenfeld2 (LR) have introduced for these systems a quantum-mechanical LR invariant, and derived a relation between the invariant eigenstates and exact solutions of the corresponding time-dependent Schrödinger equation. Using the LR invariant method, exact wavefunctions have been obtained for harmonic oscillators with time-dependent frequency,3-5 time-dependent mass and frequency,6,7 and linear driving terms.8 At present, the exact solution can be, in principle, constructed for a general one-dimensional time-dependent Hamiltonian of $`N`$ coupled quantum oscillators.9,10 The addition of a singular inverse quadratic potential does not break the exact solvability of the time-dependent harmonic oscillator, as has been shown by combining the LR invariant method with a unitary transformation in Refs. 14-17.

In this paper, we present a new vast array of quantum models with time-dependent Hamiltonians of a special scaling form. By using a couple of time-dependent unitary transformations, time evolution of these models is expressed in terms of related systems with time-independent Hamiltonians. The mapping of dynamics can be performed in any dimension, for an arbitrary number of interacting particles and for any type of the interaction potential. The exact solvability of the time-independent Hamiltonian automatically means the exact solvability of the original problem with model time-dependence; otherwise, dynamics induced by the time-dependent Hamiltonian can be studied by simpler techniques known for time-independent Hamiltonians.

We first restrict ourselves to the case of one particle situated in a $`d`$-dimensional space of points $`𝐫=(x_1,x_2,\mathrm{},x_d)`$. The proposed time-dependent Hamiltonian has the following scaling form:

$$\widehat{𝐇}(t)=\frac{f^{}(t)}{f^{}(0)}\left[\frac{f(t)}{f(0)}\frac{\widehat{𝐩}^2}{2m}+\frac{f(0)}{f(t)}V\left(\frac{f(0)}{f(t)}𝐫\right)\right],$$
(1)
where $`\widehat{𝐩}=\mathrm{i}\mathrm{}_𝐫`$ is the momentum operator, $`V`$ is an arbitrary one-body potential and the prime denotes the differentiation with respect to the argument. The normalization of parameters was chosen such that at initial time $`t=0`$ the Hamiltonian takes the standard form

$$\widehat{𝐇}(t=0)=\frac{\widehat{𝐩}^2}{2m}+V(𝐫).$$
(2)
The function $`f(t)`$ is assumed to be real continuous function of time in order to ensure the hermiticity of the Hamiltonian. To simplify the notation, without any lost of generality we shall assume that

$$f(0)=1,f^{}(0)=ϵ(ϵ\mathrm{real}\mathrm{and}\mathrm{nonzero})$$
(3)
and rewrite the Hamiltonian (1) as follows

$$\widehat{𝐇}(t)=\frac{f^{}(t)}{ϵ}\left[f(t)\frac{\widehat{𝐩}^2}{2m}+\frac{1}{f(t)}V\left(\frac{𝐫}{f(t)}\right)\right].$$
$`(1^{})`$
Two choices of the function $`f(t)`$ are of special interest. The first one,

$$f(t)=(1+2ϵt)^{1/2}$$
(4)
with $`ϵ>0`$ in order to avoid singularities in time, corresponds to the constant particle mass. The corresponding Hamiltonian (1) takes the form

$$\widehat{𝐇}(t)=\frac{\widehat{𝐩}^2}{2m}+\frac{1}{1+2ϵt}V\left(\frac{𝐫}{(1+2ϵt)^{1/2}}\right).$$
(5)
The classical version of this Hamiltonian has been introduced in Ref. 18 in connection with the study of adiabatic propagation of distributions, with $`ϵ`$ taken as the adiabatic slowness parameter. The second choice is

$$f(t)=\mathrm{exp}(ϵt),$$
(6)
and the corresponding Hamiltonian reads

$$\widehat{𝐇}(t)=\mathrm{e}^{2ϵt}\frac{\widehat{𝐩}^2}{2m}+V\left(\mathrm{e}^{ϵt}𝐫\right).$$
(7)
In one dimension and when $`V(x)=m\omega ^2x^2/2`$, the model (7) is known as the Caldirola-Kanai oscillator.19-21 Its exact quantum states are known, even in a more general case with the singular inverse square potential16, $`V(x)=m\omega ^2x^2/2+k/x^2`$.

The evolution of the system with Hamiltonian $`(1^{})`$ is governed by the time-dependent Schrödinger equation

$$\mathrm{i}\mathrm{}\frac{}{t}\psi (𝐫,t)=\widehat{𝐇}(t)\psi (𝐫,t),$$
(8)
with the initial condition for the wavefunction

$$\psi (𝐫,t=0)=\psi (𝐫).$$
(9)
In the next two paragraphs, we apply successively two time-dependent unitary transformations,

$$\{\widehat{𝐇}(t),\psi (𝐫,t)\}\{\widehat{𝐇}_1,\psi _1(𝐫,t)\}\{\widehat{𝐇}_2,\psi _2(𝐫,t)\},$$

to go from the original time-dependent Hamiltonian $`\widehat{𝐇}(t)`$ to the new time-independent Hamiltonian $`\widehat{𝐇}_2`$, via an intermediate one $`\widehat{𝐇}_1`$.

We start with the unitary transformation

$`t^{}`$ $`=`$ $`{\displaystyle \frac{1}{ϵ}}\mathrm{ln}f(t),`$ (10a)
$`x_j^{}`$ $`=`$ $`{\displaystyle \frac{x_j}{f(t)}},j=1,2,\mathrm{},d.`$ (10b)
It is easy to verify that it holds

$`{\displaystyle \frac{}{t}}`$ $`=`$ $`{\displaystyle \frac{f^{}(t)}{f(t)}}\left[{\displaystyle \frac{1}{ϵ}}{\displaystyle \frac{}{t^{}}}𝐫^{}_𝐫^{}\right],`$ (11a)
$`{\displaystyle \frac{^2}{x_j^2}}`$ $`=`$ $`{\displaystyle \frac{1}{[f(t)]^2}}{\displaystyle \frac{^2}{x_{j}^{}{}_{}{}^{2}}},j=1,2,\mathrm{},d.`$ (11b)
The transformation (10) converts the evolution equation (8) to the new one

$`\mathrm{i}\mathrm{}{\displaystyle \frac{}{t}}\psi _1(𝐫,t)`$ $`=`$ $`\widehat{𝐇}_1\psi _1(𝐫,t),`$ (12a)
$`\widehat{𝐇}_1`$ $`=`$ $`{\displaystyle \frac{\widehat{𝐩}^2}{2m}}ϵ𝐫\widehat{𝐩}+V(𝐫).`$ (12b)
The original wavefunction $`\psi `$ in (8) is expressible in terms of $`\psi _1`$ as follows

$$\psi (𝐫,t)=\psi _1(\frac{𝐫}{f(t)},\frac{1}{ϵ}\mathrm{ln}f(t)),$$
(13)
and the initial condition (9) now reads

$$\psi _1(𝐫,t=0)=\psi (𝐫).$$
(14)
The time-independent Hamiltonian $`\widehat{𝐇}_1`$ in (12) can be reexpressed as follows

$`\widehat{𝐇}_1`$ $`=`$ $`{\displaystyle \frac{1}{2m}}\left(\widehat{𝐩}mϵ𝐫\right)^2+{\displaystyle \frac{1}{2}}ϵ\left(\widehat{𝐩}\widehat{𝐫}\widehat{𝐫}\widehat{𝐩}\right)`$ (15)
$`+V(𝐫){\displaystyle \frac{1}{2}}mϵ^2r^2.`$
Since

$$\widehat{𝐩}\widehat{𝐫}=\widehat{𝐫}\widehat{𝐩}\mathrm{i}\mathrm{}d,$$
(16)
the evolution Eq. (12) can be written as

$`\mathrm{i}\mathrm{}\left({\displaystyle \frac{}{t}}+{\displaystyle \frac{dϵ}{2}}\right)\psi _1(𝐫,t)`$ $`=`$ $`{\displaystyle \frac{1}{2m}}\left(\widehat{𝐩}mϵ𝐫\right)^2\psi _1(𝐫,t)`$ (17)
$`+\left[V(𝐫){\displaystyle \frac{1}{2}}mϵ^2r^2\right]\psi _1(𝐫,t).`$
Using the unitary transformation

$$\psi _1(𝐫,t)=\mathrm{exp}\left(\frac{dϵ}{2}t+\frac{\mathrm{i}mϵ}{2\mathrm{}}r^2\right)\psi _2(𝐫,t)$$
(18)
we finally arrive at

$`\mathrm{i}\mathrm{}{\displaystyle \frac{}{t}}\psi _2(𝐫,t)`$ $`=`$ $`\widehat{𝐇}_2\psi _2(𝐫,t),`$ (19a)
$`\widehat{𝐇}_2`$ $`=`$ $`{\displaystyle \frac{\widehat{𝐩}^2}{2m}}+\left[V(𝐫){\displaystyle \frac{1}{2}}mϵ^2r^2\right].`$ (19b)
The initial condition for $`\psi _2`$ is implied by Eqs. (14) and (18) as follows

$$\psi _2(𝐫,t=0)=\psi (𝐫)\mathrm{exp}\left(\frac{\mathrm{i}mϵ}{2\mathrm{}}r^2\right).$$
(20)
The original wavefunction $`\psi `$ is obtained in terms of $`\psi _2`$ by combining relations (13) and (18),

$$\psi (𝐫,t)=\frac{1}{[f(t)]^{d/2}}\mathrm{exp}\left(\frac{\mathrm{i}mϵ}{2\mathrm{}}\left[\frac{r}{f(t)}\right]^2\right)\psi _2(\frac{𝐫}{f(t)},\frac{1}{ϵ}\mathrm{ln}f(t)).$$
(21)
Using the above mapping of dynamics, standard techniques able to deal with time-independent Hamiltonians (see e.g. Ref. 22) can be directly applied to time-dependent Hamiltonians. The mapping becomes even more appealing when the resulting time-independent problem with potential

$$V(𝐫)\frac{1}{2}mϵ^2r^2$$

is exactly solvable, because this automatically means the exact solvability of the original model time-dependence. It is clear that $`V(𝐫)`$ must be expressible as an integrable potential plus the harmonic-oscillator potential $`mϵ^2r^2/2`$. In one dimension, there exists an infinite chain of exactly solvable “reflectionless potentials”.23,24 In higher dimensions, there are standard models which admit the exact solution.22 We briefly discuss the simplest example – the isotropic $`d`$-dimensional harmonic oscillator with $`V(𝐫)=m\omega ^2r^2/2`$. For the special choices (4) and (6) of the function $`f(t)`$, the respective time-dependent Hamiltonians

$$\widehat{𝐇}(t)=\frac{\widehat{𝐩}^2}{2m}+\frac{1}{2}m\left(\frac{\omega }{1+2ϵt}\right)^2r^2$$
(22)
and

$$\widehat{𝐇}(t)=\mathrm{e}^{2ϵt}\frac{\widehat{𝐩}^2}{2m}+\frac{1}{2}\mathrm{e}^{2ϵt}m\omega ^2r^2$$
(23)
are transformed to the same time-independent Hamiltonian

$$\widehat{𝐇}_2=\frac{\widehat{𝐩}^2}{2m}+\frac{1}{2}m(\omega ^2ϵ^2)r^2.$$
(24)
The time evolution associated with this Hamiltonian can be obtained by using the method of the separation of variables. Without going into details, there are three regimes: if $`\omega >ϵ`$ one expands the wavefunction $`\psi _2`$ (and, consequently, $`\psi `$) in Hermite polynomials, if $`\omega =ϵ`$ in plane waves and if $`\omega <ϵ`$ in hypergeometric functions.

The extension of the formalism to the case of an arbitrary number of particles is straightforward and we only write down final formulae. For $`N`$ particles, the Hamiltonian $`(1^{})`$ generalizes to

$$\widehat{𝐇}(t)=\frac{f^{}(t)}{ϵ}\left[f(t)\underset{j=1}{\overset{N}{}}\frac{\widehat{𝐩}^2}{2m_j}+\frac{1}{f(t)}V\left(\left\{\frac{𝐫_j}{f(t)}\right\}\right)\right],$$
(25)
where $`f(0)=1`$ and $`f^{}(0)=ϵ`$, and the interaction function $`V`$ \[with particle coordinates uniformly scaled by $`1/f(t)`$\] involves all possible one-, two-, $`\mathrm{}`$, $`N`$-body potentials. The solution of the time-dependent Schrödinger equation

$$\mathrm{i}\mathrm{}\frac{}{t}\psi (\{𝐫_j\},t)=\widehat{𝐇}(t)\psi (\{𝐫_j\},t)$$
(26)
with the initial condition

$$\psi (\{𝐫_j\},t=0)=\psi (\{𝐫_j\})$$
(27)
reads

$$\psi (\{𝐫_j\},t)=\frac{1}{[f(t)]^{Nd/2}}\mathrm{exp}\left(\frac{\mathrm{i}ϵ}{2\mathrm{}}\underset{j=1}{\overset{N}{}}m_j\left[\frac{r_j}{f(t)}\right]^2\right)\psi _2(\left\{\frac{𝐫_j}{f(t)}\right\},\frac{1}{ϵ}\mathrm{ln}f(t)).$$
(28)
Here, the wavefunction $`\psi _2`$ evolves under the action of a time-independent Hamiltonian,

$`\mathrm{i}\mathrm{}{\displaystyle \frac{}{t}}\psi _2(\{𝐫_j\},t)`$ $`=`$ $`\widehat{𝐇}_2\psi _2(\{𝐫_j\},t),`$ (29a)
$`\widehat{𝐇}_2`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \frac{\widehat{𝐩}_j^2}{2m_j}}+\left[V(𝐫){\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \frac{1}{2}}m_jϵ^2r_j^2\right],`$ (29b)
with the initial condition

$$\psi _2(\{𝐫_j\},t=0)=\psi (\{𝐫_j\})\mathrm{exp}\left(\frac{\mathrm{i}ϵ}{2\mathrm{}}\underset{j=1}{\overset{N}{}}m_jr_j^2\right).$$
(30)
As the simplest exactly solvable example, we mention a system of $`N`$ coupled $`d`$-dimensional harmonic oscillators, which is solvable in the transformed picture with the time-independent Hamiltonian by the ordinary normal-mode technique.

The couple of time-dependent unitary transformations presented in this paper can be applied directly to the quantum Liouville equation for the density matrix, extending in this way the treatment also to the propagation of distributions.

In conclusion, we believe that the results presented will enable one to answer some problems concerning the evolution of quantum systems under a time-dependent Hamiltonian, for instance, to predict a qualitative change in the analytic structure of leading corrections to the ideal adiabaticity when going from a few degrees of freedom to a dissipative system with infinite degrees of freedom. For classical adiabatic processes in one dimension, the corresponding exact analysis has been done in Ref. 17.

Acknowledgements: I am grateful to Prof. J. K. Percus for stimulating discussions. This work was supported by Grant VEGA 2/7174/20.


text:
# Algorithm for Linear Response Functions at Finite Temperatures: Application to ESR spectrum of 𝑠=1/2 Antiferromagnet Cu benzoate

## I  Introduction

Computational physicists often face to problems of calculating the linear response functions or the real-time Green’s functions of quantum manybody systems with degree of freedom $`N10^6`$ or more. Direct diagonalization of the Hamiltonian matrix is the most naive and powerful method for modest system size $`N10^3`$. However, it becomes prohibitive for larger system because the computational time grows as $`N^3`$. Therefore efficient numerical algorithms, such as Quantum Monte Carlo methods, Lanczos Methods, and Kernel Polynomial Methods have been developed and applied to various problems.

Quantum Monte Carlo methodsSuzuki1977 ; Linden1992  can generate the Green’s functions of very large systems since QMC does not need to store the wave functions. They have been successfully used for evaluating the imaginary-time Green’s functions and, therefore, various thermodynamic quantities. However, for evaluating dynamical quantities such as AC conductivity, one has to rely on numerical analytic continuation (e.g. Maximum Entropy Method Silver1990  ) to obtain the real-time Green’s function from the imaginary-time one. This procedure is not straightforward because the statistical errors are amplified by numerical analytical continuation and the default model for the MEM must be assumed a priori.

Lanczos MethodsLanczos1950 ; Dagotto1994  have been useful techniques for evaluating dynamical responses of relatively large systems. The LM uses Lanczos recursion formula with Matrix Vector Multiplications ($`|\varphi ^{}=H|\varphi `$) to tridiagonalize the Hamiltonian matrix, which leads to a continued fraction representation of the real-time Green’s function. The drawback of LM is numerical instability of Lanczos recursion formula caused by accumulation of roundoff errors for large numbers of MVM’s. Recently, LM was extended to finite temperatures (Finite Temperature LM) by introducing random sampling over the ground and excited eigenstates Jaklic1994 . However, FTLM has a weak point that the number of excited eigenstates to be calculated increases rapidly as temperature or system size becomes large. Therefore reduction of computational costs by exploiting symmetries of the system becomes crucial.

Kernel Polynomial Method Silver1994 ; Wang1994  calculates density of states and linear response functions by using Chebyshev polynomial expansion Kosloff1986 ; Vijay2002 ; Recipe  of broadened delta functions. The Chebyshev polynomials are obtained through Chebyshev recursion formula by using MVM’s. Unlike Lanczos recursion, Chebyshev recursion avoids accumulation of numerical roundoff errors even for large numbers of MVM’s. The use of KPM for linear response functions has been limited to the ground state calculation.

Fermi-Weighted Time-Dependent Method Iitaka1997  is a combination of conventional time-dependent method Heller1993 ; Iitaka1994 ; Williams1985 , random vector representation of trace Hams2000 , and Chebyshev polynomial expansion of step function $`\theta (E_fH)`$ to extract occupied and unoccupied one-particle states of noninteracting fermi particles in the ground state. This method was successfully applied to optical absorption of silicon nano-crystallitesNomura1997 .

Our new algorithm, the Boltzmann-Weighted Time-Dependent Method, calculates linear response functions at finite temperatures by using the Boltzmann-weight operator in place of the Fermi-weight operator. The BWTDM has an advantage to FTLM that the computational time does not increase as temperature increases, because BWTDM does not calculate any eigenstates at all. Therefore BWTDM does not need use of symmetries and can easily treat disordered systems.

In this Letter, we introduce BWTDM and apply it to study Electron Spin Resonance spectrum of one-dimensional $`s=\frac{1}{2}`$ antiferromagnet Cu benzoate $`Cu(C_6H_5COO)_23H_2O`$, especially the dynamical crossover between spinon excitation and breather excitation as a function of temperature. Nonperturbative analytical calculation of this phenomena is a very difficult problem which still open. Even within Sine Gordon theory, the linear response functions have been studied only at high temperature regime (by perturbation) and at zero temperature. We compare our numerical results with preceding experimental and theoretical studies Oshima1978 ; Asano2000 ; Ogasahara2002 ; Oshikawa1999 ; Affleck1999 ; Oshikawa2002 ; Oshikawa2002b ; Oshikawa1997 ; Essler1998 ; Essler1999 ; Lou2002 .

## II The Model

The effective Hamiltonian of one-dimensional $`s=\frac{1}{2}`$ antiferromagnet Cu benzoate may be written Oshikawa1997 ; Essler1998 ; Essler1999  as

$`H`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{N_{spin}}{}}}\left[J𝐬_j𝐬_{j+1}g\mu _BH_xs_j^x+(1)^jg\mu _Bhs_j^y\right].`$ (1)
The first term is the isotropic Heisenberg antiferromagnet where $`𝐬_j`$ is spin operator at the $`j`$-th site, $`J/k_B=17.2K=1.48(meV)`$ is the exchange interaction Oshima1978 , and the sums are taken over $`N_{spin}`$ spins with periodic boundary conditions. The second is the normal Zeeman term due to the external uniform magnetic field $`H_x`$. The third is the staggered Zeeman term due to induced staggered magnetic field $`h=0.065H_x`$ originating from the Dzyaloshinskii-Moriya (DM) interaction and the staggered component of the $`g`$ tensorOshikawa2002 . In the following we assume that $`g=2.25(2.29)`$ with $`H_x`$ along $`c`$($`c^{\prime \prime }`$)-direction. Unit of magnetic field and frequency become $`J/g_c\mu _B=11.4(T)`$, $`J/g_{c^{\prime \prime }}\mu _B=11.2(T)`$, and $`J/h=359(GHz)`$, respectively.

According to the linear response theory, absorption of electromagnetic waves of frequency $`\omega `$ and polarization $`\mu `$ is expressed as

$$I^\mu (\omega )=\frac{H_R^2\omega }{2}\chi _{\mu \mu }^{^{\prime \prime }}(q=0,\omega )$$
(2)
where $`H_R`$ is the amplitude of the electromagnetic waves and $`\chi _{\mu \mu }^{^{\prime \prime }}(q,\omega )`$ is the imaginary part of the dynamical magnetic susceptibility,

$$\chi _{\mu \mu }^{^{\prime \prime }}(q,\omega )=(1e^{\beta \omega })\mathrm{Im}\underset{\eta +0}{lim}_0^{\mathrm{}}𝑑te^{i(\omega i\eta )t}g_q^\mu (t)$$
(3)
at inverse temperature $`\beta =1/(k_BT)`$. The correlation function is defined by

$$g_q^\mu (t)=\mathrm{Tr}\left[e^{\beta H}M_q^\mu e^{+iHt}M_{+q}^\mu e^{iHt}\right]/\mathrm{Tr}\left[e^{\beta H}\right]$$
(4)
where the magnetization operator is $`M_q^\mu =_{j=1}^{N_{spin}}s_j^\mu \frac{e^{iqj}}{\sqrt{N_{spin}}}.`$

## III Algorithm

The essence of BWTDM is evaluation of eq. (4) by using Boltzmann-weighted random vectors, $`|\mathrm{\Phi }_{Boltz}=e^{\beta H/2}|\mathrm{\Phi }`$,

$`g_q^\mu (t)`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Phi }_{Boltz}|M_q^\mu e^{+iHt}|M_{+q}^\mu |e^{iHt}|\mathrm{\Phi }_{Boltz}}{\mathrm{\Phi }_{Boltz}|\mathrm{\Phi }_{Boltz}}}.`$ (5)
In the above equation, the uniform random vector, $`|\mathrm{\Phi }`$, is defined by

$$|\mathrm{\Phi }\underset{n=1}{\overset{N}{}}|n\xi _n$$
(6)
where $`\{|n\}`$ is the basis set composed of eigenvectors of $`s_j^z`$’s, and $`\xi _n`$ are a set of complex random variables satisfying

$$\xi _n^{}^{}\xi _n=\delta _{n^{}n}.$$
(7)
Here $``$ stands for the statistical average. Then $`\mathrm{\Phi }|X|\mathrm{\Phi }`$ gives the trace of $`X`$,

$`\mathrm{\Phi }|X|\mathrm{\Phi }={\displaystyle \underset{n}{}}\xi _n^{}\xi _nn|X|n`$ (8)
$`+{\displaystyle \underset{nn^{}}{}}\xi _n^{}^{}\xi _nn^{}|X|n`$
$`=`$ $`{\displaystyle \underset{n}{}}n|X|n=\mathrm{tr}\left[X\right]`$
where we used eq. (7). The Boltzmann operator is evaluated by the Chebyshev polynomial expansion Kosloff1986 ; Vijay2002 ; Recipe ,

$$e^{\tau H}|\varphi =\underset{k}{}a_k(\tau )T_k(H)|\varphi $$
(9)
where $`a_k(\tau )`$ is the coefficient of the Chebyshev polynomial expansion, and $`|\varphi `$ is an arbitrary vector. The Chebyshev polynomial of k-th order $`T_k(x)`$ is calculated by the Chebyshev recursion formula,

$$T_{k+1}(H)|\varphi =2HT_k(H)|\varphi T_{k1}(H)|\varphi .$$
(10)
The time-evolution operator $`|\varphi ;t+\mathrm{\Delta }t=e^{iH\mathrm{\Delta }t}|\varphi ;t`$ can be evaluated with numerical stability by the leap frog method Iitaka1994 ,

$$|\varphi ;t+\mathrm{\Delta }t=2iH\mathrm{\Delta }t|\varphi ;t+|\varphi ;t\mathrm{\Delta }.$$
(11)
In real calculation, we first generate $`N`$ complex random numbers $`\xi _n`$ and construct $`|\mathrm{\Phi }`$ according to eq. (6). The $`|\mathrm{\Phi }_{Boltz}`$ can be computed with numerical stability by applying eq. (9) repeatedly to $`|\mathrm{\Phi }`$. Then $`\mathrm{\Phi }_{Boltz}|\mathrm{\Phi }_{Boltz}`$ gives a sample of the denominator of eq. (5). Here, we introduce another vector $`|\mathrm{\Phi }_{M_{+q}^\mu }=M_{+q}^\mu |\mathrm{\Phi }_{Boltz}`$ and calculate the time evolution of $`|\mathrm{\Phi }_{Boltz}`$ and $`|\mathrm{\Phi }_{M_{+q}^\mu }`$ by eq.(11). At each time $`t=n\mathrm{\Delta }t`$, the matrix element $`\mathrm{\Phi }_{M_{+q}^\mu };t|M_{+q}^\mu |\mathrm{\Phi }_{Boltz};t`$ gives a sample of numerator of eq. (5). After calculating the denominator and numerator for $`N_{rand}`$ random vectors, the ratio of their averages gives $`g_q^\mu (t)`$ of eq. (5). Then $`\chi _{\mu \mu }^{^{\prime \prime }}(q,\omega )`$ with a finite frequency resolution $`\eta `$ is calculated by using Gaussian filter,

$$\chi _{\mu \mu }^{^{\prime \prime }}(q,\omega )=(1e^{\beta \omega })\mathrm{Im}_0^{T_{max}}𝑑te^{i\omega t}g_q^\mu (t)\left[Ae^{\eta ^2t^2/2}\right]$$
(12)
where $`T_{max}`$ is related to $`\eta `$ by $`T_{max}\eta ^1`$. In order to avoid finite size effect of spinon excitations, $`T_{max}`$ is chosen so that $`T_{max}<N_{spin}/v_{spin}`$ where $`v_{spin}1`$ is the spin wave velocity Essler1999 . This limits the finest frequency resolution. Note that this is not the limit due to the algorithm but due to the finite size model. Much finer resolution can be used for breather modes, which are spatially localized.

Two more numerical details: First, the fluctuation in the result is inversely proportional to the square root of the number of eigenstates participating in $`|\mathrm{\Phi }_{Boltz}`$ multiplied by $`N_{rand}`$Hams2000 . Therefore systems at lower temperatures require more random vectors for the same quality of the result. Second, since the Chebyshev polynomials are defined only in the range $`[1,+1]`$, the Hamiltonian should be rescaled when eq. (9) is applied, so that all eigenvalues of $`H`$ lie in the range $`[1,+1]`$.

## IV Results and Discussions

Figure 1 shows the ESR spectra of normal polarization (Faraday configuration) $`I^z(\omega )`$ calculated with various temperatures, and Fig. 2 shows the width and peak frequency of the spinon excitation $`S`$ and the first breather excitation $`B_1`$ obtained by fitting two Gaussian peaks to the calculated spectra with the least square method. At high temperatures ($`\beta 5`$), $`S`$ is outstanding in the spectrum, and its peak shifts to higher frequency and becomes broadened as temperature decreases while $`B_1`$ becomes narrower. At low temperatures ($`\beta 5`$), $`B_1`$ prevails while its peak frequency is almost constant and its width becomes much smaller than the numerical resolution $`\eta `$. This crossover behavior between spinon excitation and first breather excitation has become computable by the invention of BWTDM, and the result is consistent with both experimental Asano2000  and field-theoretical results Oshikawa2002 . In addition, we can observe in Fig.2 the shift of $`B_1`$ at high temperatures and the narrowing of $`S`$ at low temperatures. We believe these tendencies are real physics but not numerical artifact though further experimental and theoretical studies are necessary to confirm it. Several other weak peaks appear in our results as well as in experimental results, which are supposed to be higher-order breathers and transitions between excited states. However, we are not going to analyze them here.

Figure 3 shows the induced energy gap $`E_g(H_x)`$ calculated from the $`B_1`$ frequency $`\omega `$ by solving Oshikawa1999

$$\mathrm{}\omega =\sqrt{(g\mu _BH_x)^2+\left(E_g(H_x)\right)^2}.$$
(13)
For $`H<1`$, $`E_g(H_x)`$ scales to $`H_x^{2/3}`$ as predicted by Oshikawa and AffleckOshikawa1997 . For $`1<H_x<4`$, our result seems to deviate from the $`H_x^{2/3}`$-scaling and have a dip around $`H_x^{2/3}2`$. This seems consistent with the experimental results and the density-matrix renormalization-group study Lou2002 . At such a high field, assumption of Sine Gordon theory, that the $`H_x`$ is weak and the frequency is small, is probably broken. Since the gap is calculated as a small difference of two large numbers, $`\mathrm{}\omega `$ and $`g\mu _BH_x`$, the agreement of our result with the experimental value and the $`H_x^{2/3}`$-scaling at low fields demonstrates the numerical accuracy of our algorithm.

In summary, we developed an efficient and stable algorithm for linear response functions at finite temperatures, and studied the ESR spectra of $`s=\frac{1}{2}`$ antiferromagnet Cu benzoate for a wide range of temperature and magnetic field. We reproduced experimental results of the spinon-breather crossover as a function of temperature. Temperature dependence of the width and shift of the peaks are also calculated, which are consistent with experiments and analytical theories. The calculated frequency of $`B_1`$ as a function of $`H_x`$ agrees well with both experimental and field-theoretical results at low fields, and reproduces the deviation of experimental results from the Sine Gordon theory at high fields, where the low energy assumption of the Sine Gordon theory may be broken and the choice of the compactification radius $`R`$ becomes ambiguousOshikawa2002b . The advantage of BWTDM is being applicable to finite temperatures, strong magnetic field and high frequency while its weak point is finite size effects ($`N_{spin}20`$). The computational cost of BWTDM is moderate. Calculation of one curve in Fig. 1, for example, requires approximately 30 minutes with 8 CPU’s of Fujitsu VPP5000 vector-parallel computer. We hope that this algorithm stimulates further numerical investigation of dynamical properties of quantum manybody systems at finite temperatures in general.

###### Acknowledgements.

One of the authors (TI) would like to thank Professor Seiji Miyashita for introducing the problem of ESR spectrum of Cu benzoate and Professor Masuo Suzuki for his continuous encouragement. The results presented here were computed by using supercomputers at RIKEN and NIG.


text:
# Charmonium near the deconfining transition on the lattice

## 1 Introduction

Investigation of the QCD phase transition at finite temperature and density is one of most important subjects in particle and nuclear physics. To create the quark gluon plasma phase (QGP), in which quarks and gluons are deconfined and the chiral symmetry is restored, the high energy heavy ion collision experiments have been performed and the RHIC experiment at BNL is in progress. The $`J/\psi `$ suppression has been considered as one of most important signals of formation of QGP phase. However, theoretical understanding of this phenomenon is still not sufficient to draw a conclusive scenario.

We investigate the charmonium properties at finite temperature using lattice QCD simulations which enable us to incorporate the nonperturbative effects of QCD from the first principle. In this paper, we report our results of three such approaches: (a) a nonrelativistic potential model analysis with the static quark potential measured in lattice QCD, (b) the spatial correlation between quark and antiquark , and (c) direct reconstruction of the spectral function from lattice data of correlators using the maximum entropy method and the standard $`\chi ^2`$ analysis . At present stage, all these works are done in the quenched approximation (without dynamical quark effects). To circumvent severe restriction of degrees of freedom in Euclidean temporal direction, we employ the anisotropic lattice on which temporal lattice spacing is finer than the spatial one (for setup of anisotropic lattice, see Ref. ).

The results of these analyses are summarized as follows. Although the simple potential model analysis (a) indicate that there is no bound state above $`T_c`$, the latter two analyses imply that the structure of QCD plasma is much involved and persistent collective modes may exist.

## 2 Potential model analysis

The static quark potential $`V_{Q\overline{Q}}(\stackrel{}{r})`$ is determined in lattice QCD simulation. It is extracted from the spatial correlation of the Polyakov loops, $`P(\stackrel{}{x})`$:

$$P(0)P^{}(\stackrel{}{r})C\mathrm{exp}(V_{Q\overline{Q}}(\stackrel{}{r})N_t),P(\stackrel{}{x})=\text{Tr}_{t=0}^{N_t1}U_4(\stackrel{}{x},t).$$
(1)
As displayed in the left panel of Fig. 1, the results above $`T_c`$ are well fitted to the Yukawa-type form, $`V_{Q\overline{Q}}(\stackrel{}{r})=A\mathrm{exp}(\mu |\stackrel{}{r}|)/r`$, with the screening length $`1/\mu `$. The obtained value of $`\mu `$ at each temperature is shown in the right panel of Fig. 1.

To study whether $`q`$-$`\overline{q}`$ system has a bound state or not with the obtained potential, we solve the Schrödinger equation for a spin averaged state, which is written for the radial wave function for the S-wave, $`\mathrm{\Phi }(x)`$, as

$$\left(\frac{\mathrm{\Delta }}{2m_R}+V_{Q\overline{Q}}(x)\right)\varphi (x)=E\varphi (x),\varphi (x)=x\mathrm{\Phi }(x),$$
(2)
where $`m_R=m_c/2`$ is the reduced mass and $`m_c=1.3`$ GeV is adopted. As the potential $`V_{Q\overline{Q}}(x)`$, we use the Yukawa-type form and obtain the upper bound of the screening mass, $`\mu _c`$, below which the bound state is formed (the solid line in Fig. 1(right)). Since the values of $`\mu `$ measured in the lattice simulation are well larger than $`\mu _c`$, we conclude that no bound state exists at $`T>1.05T_c`$ from the the simple potential model analysis with the static quark potential obtained in the quenched lattice QCD. This result is consistent with Ref. .

## 3 Spatial $`q`$-$`\overline{q}`$ correlation

We examine a spatial $`q`$-$`\overline{q}`$ correlation (Bethe-Salpeter amplitude) of the charmonium at finite temperature ,

$$w_\mathrm{\Gamma }(r,t)=_\stackrel{}{x}\overline{q}(\stackrel{}{x}+\stackrel{}{r},t)\mathrm{\Gamma }q(\stackrel{}{x},t)O^{}(0),O(\stackrel{}{x},t)=_\stackrel{}{y}\varphi (\stackrel{}{y})\overline{q}(\stackrel{}{x}+\stackrel{}{y},t)\mathrm{\Gamma }q(\stackrel{}{x},t),$$
(3)
where $`\mathrm{\Gamma }=\gamma _5`$ for PS and $`\gamma _i`$ for V channels. $`\varphi (\stackrel{}{y})`$ is a smearing function for which we use the result of $`w_\mathrm{\Gamma }(\stackrel{}{r},t)`$ at $`t1`$ at $`T=0`$. Figure 2 shows the result of $`q`$-$`\overline{q}`$ correlation in the vector channel normalized at the spatial origin, $`\varphi _\mathrm{\Gamma }(r,t)=w_\mathrm{\Gamma }(r,t)/w_\mathrm{\Gamma }(0,t)`$, together with the cases of free quark and antiquark. The $`t`$-dependence $`\varphi _\mathrm{\Gamma }(r,t)`$ give us a hint on the behavior of quark and antiquark in plasma phase, since if there is no bound state this spatial correlation becomes broader as $`t`$ increases, like the case of free quarks. Contrary to naive expectation, even at $`T1.5T_c`$ we observe strong spatial correlation similar to that below $`T_c`$, which is stable against $`t`$. This implies that the quark and antiquark behave differently from the free quarks, and still tend to stay together up to this temperature.

## 4 Spectral function

The spectral function contains the most direct information on the structure of the mesonic correlator, $`C(t)=_\stackrel{}{x}O(\stackrel{}{x},t)O^{}(0,0)`$. Recently the maximum entropy method (MEM) has been successfully applied to the extraction of the spectral function from lattice data at $`T=0`$ without assuming specific form . However, its applicability to problems at $`T>0`$ is not straightforward, because available physical range of $`C(t)`$ as well as number of data points are severely limited. We propose to use a standard $`\chi ^2`$ fit analysis in combination with MEM, since the former gives more quantitative result for the parameters such as mass and width than the latter if presumable form of the spectral function is given, which can be estimated with MEM. The numerical simulation is performed on lattices with spatial cutoff $`a_\sigma ^12`$ GeV and the anisotropy $`a_\sigma /a_\tau =4`$.

To verify the applicability of these methods at $`T>0`$, we require that they correctly work for the correlators at $`T=0`$ with restricted number of data points corresponding to the case at $`T>0`$. Unfortunately, MEM does not pass this test for the correlators between local operators with present level of statistics ($`O(500)`$ configurations): MEM fails to reproduce the correct feature of the spectral function at low energy even qualitatively for the region of $`C(t)`$ less than $`O(0.5\text{fm})`$. The reason may be that one need to extract the low energy structure from the correlator in a short distance, where the correlator contains the contribution from wide range of frequency of the spectral function. This result warns against physical significance of MEM analysis of such correlators at $`T>0`$.

To circumvent this problem, we apply the smearing technique which enhances the low energy part of the spectral function. For the correlators smeared with the wave function at $`T=0`$, both the procedures works well in the low energy region. Although the features of the collective mode such as the mass and width are unchanged by the smearing, the smearing can produce a mimic peak structure in the spectral function . In order to examine this possibility we also observe the correlators smeared with a narrower function (half-smeared correlators). The dependence of extracted spectral function on the smearing function is compared with that of correlators composed of free quarks.

At $`T0.9T_c`$ the results indicate that the structure of the ground state is almost the same as at $`T=0`$: almost the same mass with the width consistent with zero.

At $`T1.1T_c`$ the result of MEM still exhibits a peak structure around the same energy region as at $`T<T_c`$. Therefore we perform the same kind of the $`\chi ^2`$ fit analysis as at $`T<T_c`$. The fit to the 2-pole ansatz and to the relativistic Breit-Wigner (BW) type ansätze give inconsistent results, and the latter fits indicate that the spectral function has a peak with almost the same mass as at $`T<T_c`$ and the width of order of 200 MeV, as shown in Fig. 3. The result depends on the smearing function only slightly, in contrast to the free quark case, and hence we conclude that this is a physically significant structure.