text:
# Spin filters with Fano dots

## I Introduction

In recent decades the electric transport through quantum dots (QD) has been extensively studied both theoretically and experimentally mesoscopic1 ; mesoscopic2 ; mesoscopic3 . As the result a comprehensive picture of a big variety of underlying physical phenomena has emerged (See e.g. ABGreview ; Areview  and references therein). Confinement of electrons in small quantum dots leads to the necessity of taking into account their Coulomb repulsion. At finite temperatures the main effect is the Coulomb blockade KMMTWW ; AL , which manifests itself in the oscillations of the dot’s conductance versus the gate voltage applied to the dot. At very low temperature in the Coulomb blockade regime the dot’s conductance increases in so-called “odd” valleys, where the number of electrons on the dot remains odd. This is a signature of the mesoscopic Kondo effect KG ; GR ; NL , which has been experimentally observed goldhaber98 ; cronenwet98 ; schmid98 . In recent years, several authors have studied the possibility of spin polarization in transport through QD’s costi ; recher ; nglee88 ; ohno  and other nanoscopic configurations like rings popp , T-shaped spin filters kiselev  and dots connected to Luttinger liquid leads schmeltzer03 . In costi ; recher ; nglee88  the QD’s are attached to two leads (right and left), which we call the substitutional configuration and the transport properties are complementary to the case studied here.

In this report we combine the standard physics of the Coulomb blockade with the additional interference of the Fano type fano . The latter is essential only in the limit of a single conducting channel, which is weakly coupled to a side quantum dot goeres00 ; THCP . The transport properties of such setup has been studied in the work by Kang et al. kang01 . It was shown that the dependence of the channel conductance on the gate voltage applied to the dot is opposite to the conductance of the dot itself. Indeed the single-particle resonance locks the channel giving rise to a resonant dip of zero conductance. Therefore, at the Coulomb blockade regime the channel conductance should reveal a sequence of sharp dips as the function of the gate voltage. Similarly at very low temperature the Kondo effect should lead to the appearance of a shallow and broad valley in the conductance versus the gate voltage curve, which corresponds to the odd number of electrons on the dot.

If the spin degeneracy is lifted, say, by an applied magnetic field, and the splitting exceeds the temperature, the resonances for spin up and spin down electrons become well separated in the absence of the Kondo effect. This gives rise to the resonance spin polarized conductance due to the effect of Fano interference.

## II The model

In what follows we are dealing essentially with the single-level Anderson impurity model hewson93 . We take advantage of recently studied side dot model kang01 ; THCP  characterized by the one-dimensional Hamiltonian:

$`H`$ $`=`$ $`t{\displaystyle \underset{n,\sigma }{}}(c_{n,\sigma }^{}c_{n1,\sigma }+c_{n,\sigma }^{}c_{n+1,\sigma })+{\displaystyle \underset{\sigma }{}}\epsilon _{0,\sigma }n_\sigma `$ (1)
$`+`$ $`{\displaystyle \underset{\sigma }{}}(Vd_\sigma ^{}c_{0,\sigma }+V^{}c_{0,\sigma }^{}d_\sigma )+Un_{}n_{},`$
where $`n_\sigma =d_\sigma ^{}d_\sigma `$. The $`c`$-operators describe the propagation of electrons in the conducting channel, while the $`d`$ operators describe the occupation of the side dot attached to site zero of the channel. The dot is modeled as an Anderson impurity with a Hubbard interaction constant $`U>0`$ hewson93 . In this paper we consider $`t=V=1`$ and $`\epsilon _F=0`$ (half-filling) for all numerical computations. We also assume that the local magnetic field induces a Zeemann splitting of the dot levels

$$\epsilon _{0,}=\epsilon _0+\mathrm{\Delta }/2,\epsilon _{0,}=\epsilon _0\mathrm{\Delta }/2,$$
(2)
and the position of the resonant level $`\epsilon _0`$ can be adjusted by an external gate voltage.

The considered model falls into the class of Anderson Hamiltonians anderson61 . Let us emphasize that the conductance depends crucially on the geometry of the coupling between the dot and the channel unlike the thermodynamic properties of the model THCP , such as the Kondo temperature $`T_K`$, etc. Below we are mostly concerned with the ground state properties at zero temperature.

## III Fano resonance for $`U=0`$

Let us consider first the case $`U=0`$. Then the dimensionless conductance per spin channel $`g_\sigma `$ is linked to the transmission coefficient $`\tau `$ by the Landauer formula

$$g_\sigma =|\tau _\sigma |^2.$$
(3)
For $`U=0`$ the spin up and spin down electrons do not interact and we can suppress the index $`\sigma `$. The transmission can be computed within the one-particle picture for an electron moving at the Fermi energy $`\epsilon _F`$:

$`\epsilon _F\varphi _i`$ $`=`$ $`t(\varphi _{n1}+\varphi _{n+1})+V^{}\phi \delta _{n0},`$ (4a)
$`\epsilon _F\phi `$ $`=`$ $`\epsilon _{0,\sigma }\phi +V\varphi _0,`$ (4b)
where $`\delta _{nm}`$ is the Kronecker symbol, $`\varphi _n`$ refers to the amplitude of a single particle at site $`n`$ in the conducting channel and $`\phi `$ is the amplitude at the side dot.

The scattering problem is solved by the ansatz

$$\varphi _{n<0}=e^{iqn}+\rho e^{iqn},\varphi _{n>0}=\tau e^{iqn},$$
(5)
with the usual dispersion relation $`\epsilon =2t\mathrm{cos}q`$. Substituting Eq. (5) in Eq. (4) we readily find

$$\tau =\frac{2it\mathrm{sin}q}{2it\mathrm{sin}q+|V|^2(\epsilon _0\epsilon _F)^1},$$
(6)
hence the conductance

$$g_\sigma =|\tau _\sigma |^2=\left[1+\frac{\mathrm{\Gamma }^2}{4(\epsilon _{0,\sigma }\epsilon _F)^2}\right]^1,$$
(7)
where $`\mathrm{\Gamma }=2|V|^2/|v_F|`$ with $`v_F=d\epsilon /dq|\epsilon _F`$ is the Fermi velocity.

On the other hand the mean occupation number on the dot at zero temperature $`n_\sigma `$ is related to the scattering phase shift by the Friedel sum rule langreth66 ; hewson93

$$n_\sigma =\frac{1}{\pi }\mathrm{Im}\mathrm{ln}(\epsilon _\mathrm{F}\epsilon _{0,\sigma }+\mathrm{i}\mathrm{\Gamma }/2)$$
(8)
suggesting the relation

$$g_\sigma =\mathrm{cos}^2\pi n_\sigma .$$
(9)
This relation has a geometric origin and actually holds for arbitrary $`U`$ provided zero temperature limit. This is in contrast to Eqs. (6,7) which are applicable only for $`U=0`$.

The symmetric form of the Fano resonance given by Eq. (7) is based on the symmetry of the coupling term in the Hamiltonian (1). For nonzero magnetic field $`\mathrm{\Delta }\mathrm{\Gamma }`$ the two Fano resonances for spin up and spin down electrons are energetically separated. Therefore, the current through the channel is completely polarized at $`\epsilon _F=\epsilon _{0,}`$ and $`\epsilon _F=\epsilon _{0,}`$.

## IV Fano resonance for $`U0`$ at the degenerate level $`\mathrm{\Delta }=0`$

Let us consider the case of zero magnetic field $`\mathrm{\Delta }=0`$ but finite Coulomb interaction $`U`$. For zero coupling, $`V=0`$, the single particle excitations on the dot are located at energies $`\epsilon _0`$ and $`\epsilon _0+U`$. When the dot is weakly connected to the conducting channel (and at very low temperatures, $`T<T_K`$), a resonant peak in the density of states on the dot arises at the Fermi energy due to the Kondo effect. By using an approach identical to that of Refs. GR ; NL , we conclude that the channel conductance $`g_\sigma `$ is controlled entirely by the expectation value of the particle number $`n_\sigma `$ as given by Eq. (9). This reasoning justifies Eq. (9) for any values of $`U`$ and $`\mathrm{\Delta }`$ kang-remark . Note that in the standard setup the dot’s conductance is proportional to $`\mathrm{sin}^2(\pi n_\sigma )`$ langreth66 ; GR ; NL  due to the different geometry of the underlying system.

If the resonant level is spin degenerate $`\mathrm{\Delta }=0`$ one has $`n_{}=n_{}`$ independent of the value of $`U`$. The mean occupation of the dot $`n_\sigma `$ is a monotonous function of $`\epsilon _0`$, ranging from 1 for $`\epsilon _0\epsilon _F`$ to 0 for $`\epsilon _0\epsilon _F+U`$. According to Eq. (9) the channel conductance vanishes for $`n_\sigma =1/2`$.

In order to compute $`n_\sigma `$ for finite $`U`$ we use the numerical techniques developed in Ref. THCP . We compare our results to those obtained from the exact solution for $`n_\sigma `$ by Wiegmann and Tsvelick WT . The latter holds only for the linearized spectrum (which is not a serious restriction at half-filling) and has a compact form only in the case of zero magnetic field.

Numerically we obtain $`n_\sigma `$ by integrating the density of states at the dot over the energy. This quantity has been previously calculated using a combined method. In the first place we consider an open finite cluster of $`N=9`$ sites which includes the impurity. This cluster is diagonalized using the Lanczos algorithm lanczos . We then embed the cluster in an external reservoir of electrons, which fixes the Fermi energy of the system, attaching two semi-infinite leads to its right and left ferrari . This is done by calculating the one-particle Green’s function $`\widehat{G}`$ of the whole system within the chain approximation of a cumulant expansion cumulant  for the dressed propagators. This leads to the Dyson equation $`\widehat{G}=\widehat{G}\widehat{g}+\widehat{T}\widehat{G}`$, where $`\widehat{g}`$ is the cluster Green’s function obtained by the Lanczos method. Following Ref. ferrari , the charge fluctuations inside the cluster are taken into account by writing $`\widehat{g}`$ as a combination of $`n`$ and $`n+1`$ particles with weights $`1p`$ and $`p`$ respectively: $`\widehat{g}=(1p)\widehat{g}_n+p\widehat{g}_{n+1}`$. The total charge of the cluster and the weight $`p`$ are calculated by solving the following equations

$`Q_c`$ $`=`$ $`n(1p)+(n+1)p,`$ (10a)
$`Q_c`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}{\displaystyle _{\mathrm{}}^{\epsilon _F}}{\displaystyle \underset{i}{}}\mathrm{Im}G_{ii}(\epsilon )d\epsilon ,`$ (10b)
where $`i`$ counts the cluster sites. Once convergence is reached, the density of states is obtained from $`\widehat{G}`$. One should stress, however, that the method is reliable only if the cluster size is larger then the Kondo length.

We also perform the alternative calculation of the conductance evaluating the local density of states $`\rho _\sigma (\epsilon )`$ at site $`0`$ meir

$$g_\sigma =2\pi t\rho _\sigma (\epsilon _F).$$
(11)
We plot $`n_\sigma `$ versus the resonance level energy $`\epsilon _0`$ for $`U=0`$ and $`U=12`$. In the presence of Coulomb repulsion the mean occupancy of the dot $`n_\sigma `$ has a flat region near the value $`1/2`$ due to the Kondo effect. This corresponds to a valley of small channel conductance (see Fig. 2). Note that similar results for the limit $`U\mathrm{}`$ have been recently obtained by Bulka and Stefanski bulka01 .

The conductance found from the numerical evaluation of $`n_\sigma `$ for finite $`U`$ is in better agreement with the exact results WT  than the conductance computed from Eq. (11). Presumably, the discrepancy originates from the finite-size effects (the Kondo length is larger than the system size used in the numerical computation). The numerical calculation of $`n_\sigma `$, which involves the energy integration, appears to be more precise then that of the local density $`\rho _\sigma (\epsilon _F)`$.

## V Spin filter for $`U0`$ and $`\mathrm{\Delta }0`$

In the presence of a magnetic field the dot occupation for spin up and spin down electrons become different. At zero temperature the conductance in each spin channel can be found from Eq. (9). The expectation value $`n_\sigma `$ decreases from $`1`$ to $`0`$ with increasing the resonant energy $`\epsilon _0`$, but it takes the value $`1/2`$ at different gate voltages for $`\sigma =`$ and $`\sigma =`$ (see Fig. 3).

Consequently we may obtain a perfect spin filter as in the case of noninteracting electrons. In Fig. 4 we plot the numerical results for the conductance obtained from Eqs. (9) and (11) in a similar manner as in Fig. 2 for $`U=0`$ and $`U=12`$ with $`\mathrm{\Delta }=1.6`$.

Since the dot levels are no longer degenerate the Kondo effect is totally suppressed even at zero temperature. Indeed no plateaus formed near $`n_\sigma =1/2`$ in Fig. 3. As the result the spin-dependent resonances in the conductance have the width $`\mathrm{\Gamma }`$. The effect of large $`U`$ manifests itself mainly in the mere enhancement of the level splitting. At low temperatures $`T\mathrm{\Delta }`$ the mean occupation number $`n_{}`$ is close to unity while the occupation number for the opposite spin direction decreases to zero. One can observe, however, from Fig. 3 that there is a mixing of the levels which gives rise to the non-monotonous behavior of $`n_\sigma `$. This results in a small dip of the spin down conductance at the resonance energy for the spin up electrons and vice versa.

## VI Conclusions

We have studied numerically the spin-dependent conductance of electrons through a single ballistic channel weakly coupled to a quantum dot. The Fano interference on the dot reverses the dependence of the channel conductance on the gate voltage compared to that of the dot’s conductance. If the spin degeneracy of the dot levels is lifted, the conductance of the channel reveals Coulomb blockade dips, or anti-resonances, as a function of the gate voltage. If the level is spin degenerate the Kondo effect suppresses the channel conductance in a broad range of the gate voltages.

When the coupling between the channel and the dot is not truly local the resonances should become highly antisymmetric. In the idealized picture of the single ballistic channel transport the zero temperature conductance vanishes at the resonance values of the gate voltage due to the Fano interference. If there is no spin degeneracy, the resonances have single-particle origin (for the simple model (1) they are separated roughly by the Coulomb energy $`U+\mathrm{\Delta }`$). At the resonance the channel conductance for a certain spin direction changes abruptly from its ballistic value $`e^2/h`$ to zero, while the opposite spin conductance remains approximately at the ballistic value.

This is in contrast with previous calculations that considered a substitutional QD, where the complementary effect is found, i.e., a perfect conductance for one spin species. However, in these cases due to the local correlations and hibridization with the leads, there remains a finite conductance of the opposite spin species blurring the filtering condition. In the lateral dot case considered here we find total suppression of one spin species, thus leading to a perfect spin filter. Moreover we obtain that the energetical splitting of the Fano resonances for spin up and spin down electrons is of the order of the Coulomb interaction $`U`$ for large $`U`$ rather than of the order of the Zeeman splitting $`\mathrm{\Delta }`$ as it would be for noninteracting electrons $`U=0`$. Thus a strong electron interaction on the dot (realistic values for $`U`$ are of the order of 1meV) will allow to apply very weak magnetic fields (less than 1Tesla) and still observe a significant splitting of the Fano resonances (provided $`k_BT<\mathrm{\Delta }`$), and thus a high quality spin filter. This is because not only the spin filter is totally reflecting electrons with the desired spin $`\sigma `$, but it will be then also close to fully transparent for electrons with opposite spin. Note, however, that the coupling to the dot should preserve the phase coherence and there should be only one conducting channel.

We thank A. Aharony, O. Entin-Wohlman, V. Fleurov, L. Glazman, K. Kikoin and M. Raikh for helpful discussions. K.H and M.E.T acknowledge CONICET for support. Part of this work was done under projects Fundación Antorchas 14116-168 and 14116-212.


format: text

text:
# Large price changes on small scales

## 1 Introduction

Research in the past years has revealed that extreme price changes are not outliers, they are significantly frequent. Analyzing how markets react to such extreme events is crucial in order to understand the price formation process on the market. Economists have analyzed how markets react to large daily price changes and significant overreaction was found in the past, although markets seem to be getting more and more efficient over time . These extreme price changes may at least partially be due to news arriving on the market, although not all events can be easily attributed to major events . Previous studies have revealed that the autocorrelation of returns on the stock markets is significant for only 15-20 minutes  thus it seems reasonable that major events take place within the trading day and daily data is not appropriate to fully understand market reaction. We will thus investigate the intraday market evolution of prices, volume and the bid-ask spread after extreme price changes of 15 minutes.

Furthermore the effect of external price shocks has been analyzed on multi-agent model simulation of stock markets  and these studies show that in case there is behavioral trading on the stock market there is overreaction to external price shocks and market volatility increases sharply at the event . Thus if we are able to localize extreme intraday price changes we will probably get an insight into the process in which the agents of different behavior form the new equilibrium market price.

## 2 Defining large intraday events

The dataset used is the TAQ database of the NYSE for the years 2000-2002. The TAQ (Trades and Quotes) Database is that supplied by the NYSE: it includes all transactions and the best bid and ask price for all stocks traded on the NYSE and on the NASDAQ. We include all stocks in our sample traded on the first trading day of 2000 both on NASDAQ and on the NYSE. A minute-to-minute dataset is generated using the last transaction, and the last bid and ask prices during every minute. Since we examine intraday price formation we will only include liquid stocks in our sample. We define liquid stocks as those for which at least one transaction was filed for at least 90% of the trading minutes of the stocks included in the DJIA during the 20 pre-event trading days.

We are studying the intraday reaction to large price shocks but we have not yet defined what we mean under the term ”large price changes on small scales”. Defining intraday 15 minute events is not an easy task because volatility is higher on average at the beginning and at the end of the day. Two trivial methods are at hand:

1. Absolute filter: using this first method we look for intraday price changes bigger than a certain level of let us say 4% within 15 minutes. In this case we have to face several problems. Most of the events we find will occur during the first or last couple of minutes of the trading day because of the U-shape intraday volatility distribution of prices. These events represent the intraday trading pattern instead of extreme events. Another problem is that a 4% price jump may be an everyday event for a volatile stock while an even smaller price move may indicate a major event in case of a low volatility stock.

2. Relative filter: in case of this second method we measure the average intraday volatility as a function of trading time during the day (this means measuring the U-shape 15-minute volatility curve for each stock prior to the event) and define an event as a price move exceeding e.g. 8 times the normal volatility during that time of the day. The problem in case of this method is the following: since price moves are very small during the noon hours, the average volatility for the 24 pre-event days in these hours may be close to zero, i.e. a small price movement (a mere shift from the bid price to the ask price for example) may be denoted as an event. In this case the events cluster around the noon hours and no events are found around the beginning and the end of the trading day.

The best solution for localizing events is an intermediary one. Let us use the relative filter and absolute filter combined: this way we can eliminate the negative effect of both filters and combine their advantages. Thus an event is taken into account if and only if it passes both the relative and the absolute filter. We adjust the absolute and relative filter so as to achieve that events are found approximately evenly distributed within the trading day. In addition we omit the first 5 minutes of trading because we do not want opening effects in our average (for market reaction at the opening minutes see ).

In order to be able to observe the exact price evolution after the intraday event it is crucial to localize the events as precisely as possible. Since some price changes may be faster than others we will allow shorter price changes than 15 minutes as well. If the price change in e.g. 8 minutes already surpasses the filter level than we will assume the price change has taken place in 8 minutes. The end of the time-window in which the event takes place is regarded as the end of the price change and this is the point to which the beginning of the post-event time scale is set: thus minute 0 is exactly the end of the earliest (and of those the shortest) time window for which the price change passes the filter. This method is constructed so as to ensure that one can definitely decide by minute 0 whether an event has taken place beforehand (in the preceding 15 minutes) or not.

We will not include events in our sample for such stocks where the price tick is high compared to nominal price: thus only stocks with a nominal price over 10 USD will be studied. Another problem arises when using NASDAQ data: there are often singular transactions filed at a price outside the bid-ask spread (sometimes even 4-8% from the mean bid-ask price). Since these do not indicate an event, they are merely an outlier they are to be excluded from the sample. Thus in all cases the trigger level specified above has to be surpassed by the change in the mean of the bid and ask prices as well (not only the change in transaction price). Choosing the sample of events included in the average is a crucial step. A major event may effect the whole market, thus a dozens of stocks may experience a price shock at the same time. Including all of the events in the average would yield a sample biased to a given major event. For convenience we will only include the events in the average which are the first during a given day for a given stock. Furthermore no overlapping events will be taken into account, thus events which happen within 15 minutes (the maximum length of an extreme price change) of the previous event taken into account in the average will be omitted.

When studying intraday price changes 1 minute volatility, 1 minute trading volume (measured in USD), and the bid-ask spread are averaged besides the cumulative abnormal returns, bid and ask prices. Volatility, trading volume, and the bid-ask spread are measured in comparison to the average minute-volatility and minute trading volume of the individual stock during the same period of the trading day (intraday volume and volatility distributions are calculated using the average of the 24 pre-event trading days). This step is crucial in order to remove the intraday U-shape pattern of the above quantities from the average since their intraday variation is in the same order of magnitude as the effect itself we are looking for .

## 3 The decay of intraday shocks

Let us turn to the results: in case of the NYSE we find 197 price increases and 182 price decreases. The number of events included in the sample are 623 and 500, respectively, on the NASDAQ. We find overreaction of around 1% both during the first 10-30 minutes after price increases and decreases (see Fig. 1), which is robust and approximately of the same size on both markets. Investigating the economic significance of these price reversals is not the scope of this paper and is addressed in a separate paper by the same authors . In case we look backwards as well, i.e. at the time before the event, significant asymmetry can be observed in the pre-event behavior before price increases and drops, both seem to occur much more frequently when the price of the stock is declining. The above asymmetry holds for NYSE data as well.

Price itself is not the only important data which gives us information on market behavior after the event. Both volatility and dollar volume increase significantly (500-800%) in comparison to their pre-event daytime-adjusted value (see Fig. 2). Volatility decreases according to a power-law on both markets (on the other hand volume does not) in the case of both increases and decreases for at least 3 decades (minute 1 to minute 1000 after the event). The exponent on the NYSE is $`0.39\pm 0.01`$ for decreases and $`0.34\pm 0.01`$ for increases, and $`0.43\pm 0.01`$ and $`0.44\pm 0.01`$, respectively, on the NASDAQ. These values indicate much faster decay of volatility than measured on autocorrelation functions of daytime-adjusted volatility for individual stocks (using the exact methodology used by Liu  for the S&P 500 stock index). Autocorrelation function of volatility exhibits power law decay as well but with a power exponent usually above $`0.2`$. Exact values on our sample period of 2000-2002 are e.g. $`0.14\pm 0.02`$ for Citigroup and $`0.20\pm 0.02`$ for Nokia. This indicates that price shocks decay much faster than usual fluctuations. A probable explanation for this phenomenon is that, as mentioned already in the introduction, price shocks are likely to be caused by external events (news), thus they can be regarded as exogenous while fluctuations in general are predominantly endogenous (and this is what is measured by the autocorrelation function). We thus propose that volatility shocks due to exogenous events decay much faster than endogenous fluctuations. Similar findings have been reported on long-term volatility behaviour after several major exogenous and endogenous daily stock price crashes . We mention here that in a recent study the bid-ask spread autocorrelations on the London stock exchange have also a power law decay . It would be interesting to compare these results with our findings too, though the effect in this quantity seems market dependent.

While the increase in volatility and volume is a robust property observable on both markets, the post-event evolution of the bid-ask spread is very different: on the NYSE the bid-ask spread widens at the event exactly like volatility and exhibits power law decay with the same exponent. On the other hand the bid-ask spread on the NASDAQ reacts to the event by a widening of a mere 20% (which is negligible in comparison to the 600% jump in the bid-ask spread on the NYSE). The cause of this is the same as which causes that there is virtually no intraday variation in the bid-ask spread on the NASDAQ : the NASDAQ is a competitive dealership market, while only one market-maker operates on the NYSE for a given stock.

## 4 Conclusions

We may thus conclude that one gets an interesting insight into the price formation process when examining large price changes on small scales. The overreaction in case of both extreme price increases and decreases on both markets refer to the presence of trading on behavioral motives within the trading day. We find that not only the volatility jumps at the event (which itself is the event) but volume, and on the NYSE the bid-ask spread increase as well. We find furthermore that the volatility decays according to a power law much faster than the autocorrelation of volatility itself. This points to the fact that extreme price shocks (which are quite probable to be exogenous) decay much faster than usual (endogenous) fluctuations of volatility.


text:
# Heating and Deceleration of GRB Fireballs by Neutron Decay

## 1 Introduction

Baryonic matter ejected in a GRB explosion is partially composed of free neutrons (see (Derishev et al., 1999) and (Beloborodov, 2003a) for detailed calculations and references). This has a strong impact on the external blast wave and changes the afterglow mechanism at radii $`r\stackrel{<}{}10^{17}`$ cm (Beloborodov, 2003b). We here investigate the dynamical effect of neutrons on the early evolution of a relativistic fireball, prior to the afterglow phase.

## 2 Fireball acceleration and neutron-proton decoupling

At the initial stage of the fireball acceleration, the neutrons are collisionally coupled to the protons and accelerated by radiation pressure at the same rate. If the energy-per-baryon, $`\eta `$, is high enough, neutrons decouple before the acceleration stage is completed. This happens if

$$\eta \eta _{}\{\begin{array}{cc}400\left[\frac{L_{52}}{R_{07}\frac{1+\xi }{2}}\right]^{1/4}\hfill & \xi 1\hfill \\ 400\left[\frac{L_{52}\xi }{R_{07}\frac{1+\xi }{2}}\right]^{1/4}\hfill & \xi >1\hfill \end{array}$$
(1)
where $`\xi `$ is the ratio of neutron and proton densities, $`R_0=10^7\mathrm{cm}R_{07}`$ is the initial size of the fireball, and $`L=10^{52}\mathrm{erg}L_{52}`$ is the equivalent isotropic luminosity of the explosion. Beyond the decoupling radius,

$$R_{np}R_0\eta _{}\left(\frac{\eta _{}}{\eta }\right)^{1/3},$$
(2)
neutrons coast with a constant Lorentz factor

$$\mathrm{\Gamma }_n\frac{R_{np}}{R_0}\eta _{}\left(\frac{\eta _{}}{\eta }\right)^{1/3}.$$
(3)
The optically thick proton fireball continues to accelerate as $`\mathrm{\Gamma }_pr/R_0`$ until the internal energy equals the proton rest-mass energy and acceleration is no longer efficient. If the fireball remains optically thick during the whole acceleration stage, its Lorentz factor saturates at

$$\mathrm{\Gamma }_p\{\begin{array}{cc}\eta \hfill & \eta <\eta _{},\hfill \\ \eta \left(1+\xi \xi \left(\frac{\eta _{}}{\eta }\right)^{4/3}\right)\hfill & \eta >\eta _{},\hfill \end{array}$$
(4)
and the acceleration stage ends at radius

$$R_s\mathrm{\Gamma }_pR_0.$$
(5)
## 3 Adiabatic cooling

We first consider the thermal history of the fireball without neutrons. Electrons, protons, and radiation maintain a common temperature $`T`$ in the early dense fireball via Coulomb collisions and Compton scattering, and $`T`$ decreases adiabatically to very low values. During such adiabatic expansion, the photon-to-baryon ratio $`\varphi =n_\gamma /n_b=const10^5`$, and radiation completely dominates the internal energy and pressure. Therefore, $`T`$ decreases according to the radiation adiabatic law with index $`\widehat{\gamma }=\frac{4}{3}`$. This continues until the fireball becomes transparent to radiation at the photosphere radius $`R_\tau `$. Assuming $`R_\tau >R_s`$, we have

$$R_\tau =\frac{L\sigma _T}{4\pi m_pc^3\mathrm{\Gamma }_p^3},$$
(6)
here $`\sigma _T`$ is Thomson cross section. After becoming transparent the plasma is still tracking the temperature of the (freely streaming) photons, which is constant in the plasma frame if the fireball coasts with a constant Lorentz factor. The electrons decouple thermally from radiation only when the Compton timescale ($`t_C=\mathrm{\Gamma }_p(3m_ec/4U_{rad}\sigma _T)`$ where $`U_{rad}`$ is the radiation energy density) exceeds the expansion timescale $`R/c`$. This happens at radius

$$R_{e\gamma }=\frac{\sigma _TL}{3\pi m_ec^3\mathrm{\Gamma }_p^{7/3}}\left(\frac{R_0}{R_\tau }\right)^{2/3}.$$
(7)
The protons decouple from the electrons when the Coulomb timescale $`t_{ep}17T^{3/2}n_e^1\mathrm{\Gamma }_p`$ s exceeds $`R/c`$. The corresponding radius of e-p decoupling is

$$R_{ep}\left(68\pi m_pc^4\right)^1\frac{L}{\mathrm{\Gamma }_p^3T_s^{3/2}}\left(\frac{R_\tau }{R_s}\right),$$
(8)
where $`T_s=T_0/\mathrm{\Gamma }_p`$ is the temperature at the saturation radius, and $`T_0`$ is the initial temperature of the fireball. The protons decouple thermally from the electrons before the electrons decouple from radiation if $`R_{ep}<R_{e\gamma }`$. One can show that $`R_{ep}/R_{e\gamma }4\times 10^2T_{s,7}^{3/2}(R_\tau /R_s)^{5/3}`$ (with $`T_s=T_{s7}\times 10^7`$ K). At radii $`R>\mathrm{min}(R_{e\gamma },R_{ep})`$, the protons are decoupled from radiation and cool adiabatically with index $`\widehat{\gamma }=\frac{5}{3}`$.

Thus, in the absence of neutrons, the fireball cools down adiabatically as $`Tr^1`$ during the acceleration stage and as $`Tr^{2/3}`$ during the subsequent coasting stage ($`\mathrm{\Gamma }_p\eta `$) up to the transparency radius. After thermal decoupling from radiation the protons cool as $`T_pr^{4/3}`$. The thermal history of a pure proton fireball is shown by dot-dashed curve in Fig. 1. For most of the early evolution, the fireball is a cold coasting outflow. The presence of neutrons and their $`\beta `$-decay change this picture.

## 4 Heating by proton–neutron collisions

Protons in the fireball are continuously heated via proton-neutron collisions (Pruet et al., 2001). The collisional heating reaches its peak at $`RR_{np}`$: at $`R<R_{np}`$ the n-p collisions are frequent but the relative velocity of the neutron and proton components is small ($`\beta _{np}\frac{t_{np}}{R/c}<1`$) and at $`R>R_{np}`$ the relative velocity is relativistic $`\beta _{np}1`$ however the collisions are rare. The collisional heating, thus, peaks at $`R_{np}`$ where the relative velocity becomes comparable to the speed of light and the rate of collisions still allows an efficient transfer of energy on a dynamical timescale. There are two sinks of heat gained by protons via n-p collisions: adiabatic cooling and Coulomb scattering off electrons. The competition between heating and cooling shapes the first peak in the proton temperature profile (Fig. 1). At $`R=R_{ep}`$ the energy transfer from the protons to the electrons becomes inefficient on the expansion timescale $`R/c`$. Then most of the heat gain by protons remains stored in the proton component, not given to the electrons and radiation.

## 5 Fireball heating and deceleration by decayed neutrons

In a fireball with $`\eta >\eta _{}`$, the neutrons have a smaller momentum than the protons after $`R_{np}`$. Their decay products $`e^{}`$ and $`p`$ exchange momentum with the ion fireball and heat it up. The heating of an ion medium by $`\beta `$-decay of neutrons moving with respect to the medium was discussed in . The decay particles exchange momentum with the medium through the two-stream instability or because of gyration in a transverse magnetic field frozen into the medium. Since the neutrons decay gradually at all radii, this momentum exchange and heating take place continuously from the moment of decoupling.

The heating rate of the fireball due to $`\beta `$-decay is

$$\frac{dE_h}{dr}=(\mathrm{\Gamma }_{rel}1)\frac{dM_{np}c^2}{dr},$$
(9)
where $`\mathrm{\Gamma }_{rel}`$ is the neutron Lorentz factor in the fireball rest frame,

$$\frac{dM_{np}}{dr}=\frac{M_n}{R_\beta },$$
(10)
is the decayed neutron mass per unit radius, and

$$R_\beta \mathrm{\hspace{0.17em}0.8}\times 10^{16}\mathrm{cm}\left(\frac{\mathrm{\Gamma }_n}{300}\right)$$
(11)
is the mean radius of decay. The heating by $`\beta `$-decay begins at $`R_{\mathrm{np}}`$ and significantly changes the thermal history of the fireball. It shapes the second peak in proton temperature (see Fig. 1).

As long as proton thermal energy is much below $`m_pc^2`$, adiabatic cooling is small compared to heating and the proton temperature increases linearly with radius,

$$T_pr.$$

At a radius $`R_{dis}`$ the dissipated bulk kinetic energy equals the total energy of the fireball and $`\mathrm{\Gamma }_p`$ begins to decrease. Between $`R_{dis}`$ and $`R_\beta `$ the $`\beta `$-decay heating balances the adiabatic cooling and

$$\mathrm{\Gamma }_p\left(\frac{R_\beta }{r}\right)^{1/2}.$$

When the fireball reaches $`R_\beta `$ the neutron component is exponentially depleted and the dissipation process switches off. The adiabatic cooling then leads to a power-law decrease in temperature,

$$T_pr^{4/3}.$$

## 6 Discussion

We have shown here that the presence of a neutron component greatly affects the early dynamics of GRB fireballs with $`\eta >\eta _{}400`$. Contrary to the pure proton model, the fireball with a neutron component heats up significantly at $`R\stackrel{>}{}R_{np}`$. In a popular GRB scenario, the bulk energy is partially converted into the observed $`\gamma `$rays through shock waves inside the fireball at $`r\stackrel{>}{}10^{12}`$ cm (Rees and Mészáros, 1994). The $`\beta `$-decay process described here should affect the internal shocks. It decelerates the fastest portions of the inhomogeneous fireball and reduces the contrast of Lorentz factors, which should significantly reduce the dissipation efficiency of internal shocks. A detailed study of these effects is in preparation.


text:
# Nearly-logarithmic decay of correlations in glass-forming liquids

(January 23, 2004)

## Abstract

Nearly-logarithmic decay of correlations, which was observed for several supercooled liquids in optical-Kerr-effect experiments \[G. Hinze et al. Phys. Rev. Lett. 84, 2437 (2000), H. Cang et al. Phys. Rev. Lett. 90, 197401 (2003)\], is explained within the mode-coupling theory for ideal glass transitions as manifestation of the $`\beta `$-peak phenomenon. A schematic model, which describes the dynamics by only two correlators, one referring to density fluctuations and the other to the reorientational fluctuations of the molecules, yields for strong rotation-translation coupling response functions in agreement with those measured for benzophenone and Salol for the time interval extending from 2 picoseconds to about 20 and 200 nanoseconds, respectively.

Optical-Kerr-effect (OKE) spectroscopy is a powerful technique for the study of the dynamics of supercooled liquids Torre et al. (1998). The experiment provides a response function $`\chi _A(t)`$ for times $`t`$ exceeding a fraction of a picosecond. The function $`\chi _A(t)`$ is proportional to the negative time derivative of a correlator $`\varphi _A(t)=A^{}(t)A/|A|^2`$. Here, $``$ denotes canonical averaging, and the probing variable $`A`$ is the anisotropic part of the dielectric function. The instrumentation was improved recently by application of heterodyne detection. As a result, it was possible to measure $`\chi _A(t)`$ for times up to 500 nanoseconds, i.e., glassy dynamics was documented for the enormous time interval starting at the end of the transient and extending over more than five orders of magnitude Hinze et al. (2000). It was shown that the evolution of the glassy dynamics of Salol upon decreasing the temperature $`T`$ can be interpreted by the universal formulas derived within the mode-coupling theory of ideal glass transitions (MCT) Torre et al. (1998); Hinze et al. (2000). The fit values for the various parameters have been found to be consistent with those obtained by other light-scattering techniques Li et al. (1992); Yang and Nelson (1995). However, the universal formulas do not describe the data for all times outside the transient regime. Rather, if $`T`$ decreases, there opens a time interval larger than two orders of magnitude, which precedes the interval of validity of the universal formulas. Here, the response for Salol follows closely a $`1/t`$-law, i.e., the correlator exhibits nearly a logarithmic decay: $`\varphi _A(t)\mathrm{ln}(t/\tau )`$ Hinze et al. (2000). This intriguing feature was observed also for some other van-der-Waals liquids, like benzophenone (BZP). This holds with the reservation that the heretofore unknown relaxation process can be described more adequately by a power law, $`\varphi _A(t)ft^b^{}`$, albeit with a rather small exponent $`b^{}`$ Cang et al. (2003). Gaining an understanding of the indicated findings is a challenge to all theories aiming to unlock a comprehensive description of liquids. In this Letter, it will be shown that the measured complex relaxation scenario Hinze et al. (2000); Cang et al. (2003) is a generic, though not universal, implication of the standard MCT.

The MCT was proposed as a mathematical model for glassy dynamics Bengtzelius et al. (1984), whose fascinating features are obtained as implication of bifurcation points of nonlinear equations of motion derived for the density fluctuations. The basic bifurcation is a fold singularity describing a transition from ergodic to non-ergodic behavior if the temperature decreases through a critical value $`T_c`$. Using $`ϵ=(T_cT)/T_c`$ as a small parameter, the long-time behavior of the correlators can be evaluated by asymptotic expansions Franosch et al. (1997). The leading-order formulas establish the universal features of the dynamics. Comparisons of experimental results with these formulas Götze (1999) and tests by molecular-dynamics simulation Götze (1999); Kob (2003) have shown that MCT is a serious candidate for an explanation of glassy dynamics. Nearly-logarithmic decay of correlations was predicted within MCT for states close to bifurcations of the cusp type as can be inferred from Ref. Götze and Sperl (2002) and the papers cited there. It was pointed out Hinze et al. (2000); Cang et al. (2003) that the OKE data might be fitted by the universal results for relaxation near cusp bifurcations. However, one has to vary at least two control parameters in order to approach such singularity. A fit of the OKE data, which only depend on the single control parameter $`T`$, would require implausible fine tuning of the coupling coefficients entering the MCT equations. The explanation proposed below is based on the existence of the simple liquid-glass-transition singularity. Our results are valid for response functions of variables $`A`$ that couple strongly to density fluctuations so that a so-called $`\beta `$-peak can occur in the susceptibility spectra Buchalla et al. (1988); Götze and Sjögren (1989).

Let us consider an MCT model describing schematically the density fluctuations by a single correlator $`\varphi (t)`$. It obeys the Zwanzig-Mori equation of motion

$$_t^2\varphi (t)+\nu _t\varphi (t)+\mathrm{\Omega }^2\varphi (t)+\mathrm{\Omega }^2_0^tm(tt^{})_t^{}\varphi (t^{})dt^{}=0,$$
(1)
where $`\mathrm{\Omega }`$ and $`\nu `$ are frequencies parameterizing the short-time asymptote, $`\varphi (t)=11/2(\mathrm{\Omega }t)^2+1/6\nu \mathrm{\Omega }^2t^3+𝒪(t^4)`$. The kernel $`m(t)`$, which represents the interactions of the density fluctuations, is modeled as

$$m(t)=v_1\varphi (t)+v_2\varphi ^2(t).$$
(2)
Here, $`v_10`$ and $`v_20`$ are the coupling coefficients. The state of the system is specified by a point $`𝐕=(v_1,v_2)`$ in the $`v_1`$-$`v_2`$ plane. For small $`𝐕`$, the model describes liquid states where $`\varphi (t\mathrm{})=0`$. For large $`𝐕`$, one gets a nontrivial long-time limit $`f=\varphi (t\mathrm{})`$, $`0<f<1`$. Parameter $`f`$ quantifies the arrest of density fluctuations in the non-ergodic glass state. There are lines of liquid-glass transition points $`𝐕^c=(v_1^c,v_2^c)`$, each characterized by some number $`\lambda ,\mathrm{\hspace{0.17em}1}/2\lambda <1`$. This number determines the critical exponent $`a`$, $`0<a<1/2`$, and the von Schweidler exponent $`b`$, $`0<b1`$, by the relation $`\mathrm{\Gamma }(1a)^2/\mathrm{\Gamma }(12a)=\lambda =\mathrm{\Gamma }(1+b)^2/\mathrm{\Gamma }(1+2b)`$. We consider the line, which can be parameterized by $`v_1^c=(2\lambda 1)/\lambda ^2,v_2^c=1/\lambda ^2,\mathrm{\hspace{0.33em}0.5}\lambda <1.`$ Crossing the line, the long-time limit of $`\varphi (t)`$ jumps from zero to the critical value $`f^c=1\lambda `$. The distance of the state $`𝐕`$ from $`𝐕^c`$ is specified by the separation parameter $`\sigma =[(v_1v_1^c)+(v_2v_2^c)f^c]f^c(1f^c),`$ which is negative for liquid states and positive for glass states. This model is the simplest one reproducing all values for the anomalous exponents, which can occur in the general theory Götze (1984). The parameter $`f^c`$ appears as plateau value of the $`\varphi (t)`$-versus-$`\mathrm{log}t`$ curves for $`𝐕`$ near $`𝐕^c`$ or as relative strength of the $`\alpha `$-peak of the loss spectrum Franosch et al. (1997).

The probing-variable correlator $`\varphi _A(t)`$ obeys the same general equation of motion and initial decay as $`\varphi (t)`$, but $`\mathrm{\Omega }`$, $`\nu `$ and $`m(t)`$ have to be replaced by the corresponding quantities $`\mathrm{\Omega }_A`$, $`\nu _A`$ and $`m_A(t)`$. The model for the kernel is

$$m_A(t)=v_A\varphi (t)\varphi _A(t),$$
(3)
where $`v_A>0`$ quantifies the coupling of the probing variable to the density fluctuations. If $`v_Af^c>1`$, the long-time limit of $`\varphi _A(t)`$ jumps from zero to the critical value $`f_A^c=11/(v_Af^c)`$ as $`𝐕`$ crosses the transition line. Originally, this model was motivated for the density fluctuations of a tagged particle Sjögren (1986). Also, the MCT equations for reorientational dynamics of a linear molecule suggest an expression like Eq. (3) Chong and Götze (2002). Within the microscopic theory, the coupling coefficients depend on the equilibrium structure functions, which depend smoothly on $`T`$. Therefore, $`v_1`$, $`v_2`$, and $`v_A`$ are smooth functions of the temperature. The specified model was applied repeatedly for the description of experiments, as can be inferred from Ref. Götze and Voigtmann (2000) and the papers cited there.

The long-time decay of the correlators at the transition point is given by $`[\varphi (t)f^c]/h=[\varphi _A(t)f_A^c]/h_A=(t_0/t)^a+𝒪(t^{2a})`$. The time $`t_0`$ depends on the transient dynamics. It is determined by matching the numerical solution of $`\varphi (t)`$ for large times to the asymptotic formula. The amplitudes are $`h=\lambda ,h_A=\lambda /(v_Af^{c\mathrm{\hspace{0.17em}2}})`$. In the limit of small $`|\sigma |`$, there appears a large time interval where $`|\varphi (t)f^c|`$ is small. Here, one gets in leading-order $`[\varphi (t)f^c]/h=[\varphi _A(t)f_A^c]/h_A=c_\sigma g(t/t_\sigma )`$. This is a scaling law with a correlation scale $`c_\sigma =\sqrt{|\sigma |}`$ and a time scale $`t_\sigma =t_0/|\sigma |^{1/2a}`$. The function $`g(\widehat{t})`$ is determined solely by $`\lambda `$. It exhibits the critical power law for small rescaled times $`\widehat{t}`$, $`g(\widehat{t}1)=1/\widehat{t}^a`$, and von Schweidler’s power law for large rescaled times in the liquid, $`g(\widehat{t}1)=B\widehat{t}^b`$ Franosch et al. (1997). The scaling law implies one for the response functions

$$\chi (t)/h=\chi _A(t)/h_A=s_\sigma k(t/t_\sigma ),$$
(4)
where $`s_\sigma =c_\sigma /t_\sigma `$. The master function $`k(\widehat{t})=\widehat{g}/\widehat{t}`$ interpolates between the critical power law, $`k(\widehat{t}1)=a/\widehat{t}^{(1+a)},`$ and, in the liquid state, the von Schweidler power law, $`k(\widehat{t}1)=bB/\widehat{t}^{(1b)}.`$

The VH-light-scattering spectra of Salol Li et al. (1992) or orthoterphenyl Cummins et al. (1997) show that $`f_A^c`$ is about $`0.9`$. This means there is strong arrest of the reorientational motion at the ideal glass-transition points in these van-der-Waals systems. Integrating $`\chi _A(t)`$ over time, one can determine $`\varphi _A(t)`$ and read off the plateau. We find $`f_A^c=0.90\pm 0.05`$ from the $`T=260\text{K}`$-data for BZP Cang et al. (2003) and $`f_A^c=0.93\pm 0.03`$ from those for Salol at $`T=257\text{K}`$ Hinze et al. (2000), corroborating the preceding conclusion. Therefore, we focus on $`\chi _A(t)`$ for large translation-rotation coupling $`v_A`$.

Figure 1 reproduces the OKE response functions for BZP Cang et al. (2003). A dashed straight line of slope $`x=0.80`$ shows that the data for $`T=251\text{K}`$ exhibit power-law decay $`\chi (t)1/t^x`$ for 2 ps$`t20`$ ps. There is a von Schweidler-law-like variation: $`\varphi (t)\text{const.}t^b^{},b^{}=1x=0.20`$. In the double-logarithmic representation, the scaling-law result, Eq. (4), appears as interpolation between a straight line of slope $`(1+a)>1`$, and one of slope $`(1b)<1`$. Changing the scales $`s_\sigma h_A`$ and $`t_\sigma `$ is equivalent to a translation of the curve. The dashed line marked sc shows an example calculated for $`\lambda =0.70`$ ($`a=0.33,b=0.64`$) shifted to match the data for $`T=251\text{K}`$. It provides a proper description of the observations for $`60\text{ps}<t<6\text{ns}`$. However, since $`1+a>1>x`$, the scaling law cannot account for the data for $`t<60`$ ps. This conclusion Cang et al. (2003) is not altered by choosing $`\lambda `$ or $`t_\sigma `$, $`s_\sigma `$ differently. The dots are the theoretical results for $`\chi _A(t)`$ scaled by amplitude factors $`N`$ and scales $`\mathrm{\Omega }`$, which vary somewhat with $`T`$. Notice that the calculated functions exhibit transient oscillations for times below and up to 2 ps, as do the measured curves. This shows that the parameters $`\mathrm{\Omega },\nu ,\mathrm{\Omega }_A,\nu _A`$ are chosen reasonably. The strong variation of $`\chi _A(t)`$ with changes of $`T`$ is due to the changes of the coupling coefficients as documented in the left inset. The state $`𝐕=(v_1,v_2)`$ shifts towards the transition line if $`T`$ decreases. An extrapolation to a transition point for $`\lambda =0.70`$, which is near the cross, is consistent with the fit values for $`v_1`$ and $`v_2`$. Parameter $`\sigma `$ exhibits an almost linear temperature dependence, and extrapolation to $`\sigma =0`$ suggests $`T_c235\text{K}`$. Uncertainty estimates for $`\lambda `$ and $`T_c`$ cannot be given yet, since we did not study in detail the possibility for data fits by other choices for $`v_1`$, $`v_2`$, and $`v_A`$.

Figure 2 shows the OKE response for Salol for the highest and lowest temperature measured and three temperatures in between Hinze et al. (2000). The data are normalized to unity for $`t=0.01`$ ps. A dashed straight line of slope $`x=1.15`$ shows that the data for $`T=247\text{K}`$ can be described by $`\chi (t)1/t^x`$ for the interval 2 ps $`<t<`$ 100 ps. Equivalently, the correlator shows a decay similar to the critical one, $`\varphi (t)1/t^a^{}`$ with $`a^{}=x1=0.15`$. For the $`T=257`$ K data, a similar behavior is found with $`x0`$, i.e., $`\varphi (t)`$ exhibits logarithmic decay Hinze et al. (2000). The MCT fits have been evaluated for temperature-independent parameters $`\mathrm{\Omega },\mathrm{\Omega }_A,\nu ,\nu _A`$. The normalization constants $`N`$ exhibit some $`T`$ dependence. The path $`𝐕=(v_1,v_2)`$ extrapolates to a critical point for $`\lambda =0.73`$ ($`a=0.31,b=0.59`$). The $`\sigma `$-versus-$`T`$ results suggest $`T_c245\text{K}`$. The numbers for $`\lambda `$ and $`T_c`$ are consistent with those obtained previously Torre et al. (1998); Hinze et al. (2000); Li et al. (1992); Yang and Nelson (1995). For the identification of the ideal transition point it would be helpful to have data available also for $`T`$ below the estimated $`T_c`$. The variation of $`v_A`$ closely follows an Arrhenius law: $`v_A\mathrm{exp}(T_0/T)`$, where the fit value $`T_0=1000\text{K}`$ is not unreasonable for a van-der-Waals liquid.

Figure 3 shows the two response functions, $`\chi (t)=\varphi (t)/t`$ and $`\chi _A(t)=\varphi _A(t)/t`$, for three liquid states near the transition point with $`\lambda =0.70`$ and the corresponding scaling-law results. There is a second scaling law, the superposition principle for the $`\alpha `$ process, dealing with the dynamics for times within the von Schweidler-law regime and longer: $`\chi (t)=\stackrel{~}{\chi }(t/t_\sigma ^{})/t_\sigma ^{},\chi _A(t)=\stackrel{~}{\chi }_A(t/t_\sigma ^{})/t_\sigma ^{}`$, $`t_\sigma ^{}=t_0/|\sigma |^\gamma ,\gamma =1/(2a)+1/(2b)`$. The insets demonstrate this result. The solution at $`𝐕^c`$ including the leading correction reads $`[\varphi _A(t)f_A^c]=h_A(t_0/t)^a[1+K_A(t_0/t)^a]`$ and a corresponding formula with index $`A`$ dropped. The number $`K_A=1.5`$ is evaluated from the coupling constants $`v_1`$, $`v_2`$, $`v_A`$; it is negative and large since $`f_A^c`$ is high Franosch et al. (1997). The result is shown as dotted line for $`\chi _A`$. The corresponding correction for $`\chi (t)`$, $`K=0.020`$, is so small that it cannot be made visible in Fig. 3. This exemplifies that the range of validity of the universal formulas can be quite different for different functions. The scaling-law formulas explain quantitatively the evolution of the glassy dynamics as exhibited by $`\chi (t)`$ and this for the complete regime outside the transient, $`\mathrm{\Omega }t>2`$. The state $`n=2`$ is marked by a cross in the left inset of Fig. 1. The states of interest for the data interpretation have a much larger separation from the critical point. One finds here, as in previous work Götze and Voigtmann (2000); Franosch et al. (1997), that the scaling-law description can explain all the qualitative features of $`\chi (t)`$, but that a quantitative description does not work for $`|\sigma |>0.05`$. The same conclusions are valid for $`\chi _A(t)`$, albeit only for times with $`t\mathrm{\Omega }3000`$. The size of the correction term $`K_A`$ is so large, that there appears a three-decade time interval outside the transient where the universal formulas cannot account for the MCT solution. Within the interval $`2\mathrm{log}_{10}(\mathrm{\Omega }t)3.5`$, the leading-order-correction formula describes the results. But for $`0.5\mathrm{log}_{10}(\mathrm{\Omega }t)2`$, even this formula is insufficient. Within the interval $`0.5\mathrm{log}_{10}(\mathrm{\Omega }t)3.5`$, nearly-logarithmic decay is exhibited by $`\varphi _A(t)`$.

The equation of motion for $`\varphi _A(t)`$ is equivalent to $`\chi _A(\omega )=\mathrm{\Omega }_A^2/[\omega (\omega +i\nu _A)\mathrm{\Omega }_A^2+\mathrm{\Omega }_A\omega m_A(\omega )]`$, where $`\chi _A(\omega )`$ and $`m_A(\omega )`$ are the Fourier transforms of $`\chi _A(t)`$ and $`m_A(t)`$, respectively, for frequency $`\omega `$. For low-frequency phenomena, $`|\omega (\omega +i\nu _A)|\mathrm{\Omega }_A^2`$, this simplifies to $`\chi _A(\omega )=1/[1\omega m_A(\omega )]`$. For $`t/t_01`$ and $`t/t_\sigma 1`$, one gets up to next-to-leading order: $`m_A(t)=f_m+h_m(t_0/t)^a[1+K_m(t_0/t)^a]`$, where $`f_m=v_Af^cf_A^c`$, $`h_m=v_A(f^ch_A+f_A^ch)`$, $`K_mh_m=v_A[f^ch_AK_A+f_A^chK+hh_A]`$. There is a trend that the negative contribution of the first term in the bracket is canceled by the positive one due to the last term. Hence, because the correction term $`K_A`$ for $`\varphi _A(t)`$ is large, the correction term $`K_m=0.16`$ for $`m_A(t)`$ is small. Thus, the critical law is a good approximation for $`m_A(t)`$. This leads to the Cole-Cole formula for the $`\beta `$-process: $`\chi _A^{\text{cc}}(\omega )=\chi _0^{\text{cc}}/[1+(i\omega /\omega _\beta )^a]`$, $`\chi _0^{\text{cc}}=1/[1+f^cf_A^cv_A]`$, $`\omega _\beta =[(1+f_m)/h_m\mathrm{\Gamma }(1a)]^{1/a}/t_0`$. The limit $`\omega /\omega _\beta 1`$ reproduces the universal critical decay: $`\chi (t)1/t^{1+a}`$. The limit $`\omega /\omega _\beta 1`$ leads to a von Schweidler-law-like decay: $`\chi (t)1/t^{1a}`$. It depends on the relative size of $`t_0^1`$, $`\omega _\beta `$ and $`t_\sigma ^1`$, which part of the $`\beta `$-peak spectrum dominates the nearly-logarithmic decay, $`\chi (t)1/t^x`$, $`1ax1+a`$. The dotted line in Fig. 3 marked cc exhibits $`\chi _A^{\text{cc}}(t)`$ for our model. The preceding reasoning for using the leading asymptotic expansion for $`m_A(t)`$ rather than for $`\varphi _A(t)`$ is valid also for the microscopic version of MCT. The schematic model used here is merely the simplest example illustrating our derivation of the $`\beta `$-peak phenomenon.

We have shown that the evolution of glassy dynamics as measured by OKE spectroscopy for two liquids can be described reasonably by the solutions of a standard schematic MCT model, and this for time intervals extending up to five orders of magnitude. The long-time part of the response functions can be explained qualitatively by the known scaling laws reflecting leading-order asymptotic solutions of the MCT equations. The leading-order correction formulas explain that for strong probing-variable-environment coupling there appears a large time interval adjacent to the one for the transient motion where the scaling laws cannot explain the data. The nearly-logarithmic decay within this interval is explained as $`\beta `$-peak dynamics and is fitted perfectly by the solutions of the schematic model. Our theory implies the prediction of nearly-logarithmic decay of reorientational correlations for liquids of molecules with large elongation. Molecular-dynamics simulations can test this and provide the structure factors for a microscopic calculation within MCT. Such work could substantiate or falsify the preceding explanation. It remains to be explained why the $`\beta `$-peak phenomenon has not been detected in earlier studies.

We thank M.D. Fayer and G. Hinze for providing their data sets. We thank them as well as H.Z. Cummins, V.N. Novikov, R. Schilling and Th. Voigtmann for stimulating discussions. Our work was supported by the DFG Grant No. Go154/13-2.


text:
# Long range Ising model for credit risk modeling in homogeneous portfolios

## Abstract

Within the framework of maximum entropy principle we show that the finite-size long-range Ising model is the adequate model for the description of homogeneous credit portfolios and the computation of credit risk when default correlations between the borrowers are included. The exact analysis of the model suggest that when the correlation increases a first-order-like transition may occur inducing a sudden risk increase. Such a feature is not reproduced by the standard models used in credit risk modeling.

preprint: ECM

The modelization and control of risk is a subject of extreme importance in a rich variety of disciplines ranging from geology to medicine. To know which is the probability distribution of rare events that may cause very important losses is the key point to be solved and to which a lot of efforts have been devoted for many years. Within the framework of finance, the most studied case is that of market risk (the risk associated to uncertainties of capital markets). But, by far, the risk that most concern practitioners (both due to its negative impact and its modeling difficulties) is credit risk T.R.Bielecki and Rutkowski (2002): a bank or financial institution give credits to $`N`$ borrowers which may fail to return it. Such a group of borrowers constitutes a so called credit portfolio. It is very important to model which is the probability $`p(L)`$ that $`L`$ ($`0LN`$) of such borrowers will fail in returning the credit. The knowledge of $`p(L)`$ is of extreme importance to decide the total amount of capital and reserves that the bank should hold to prevent a failure and for economic and regulatory reasons. Also $`p(L)`$ is the main driver in the valuation and hedging of credit derivatives.

Banks have been using many alternative models. Most of these, as a first step, consider the case of homogeneous portfolios for which the properties of the borrowers and their credits are supposed to be the same. In this paper we will also restrict to this homogeneous case. In the financial literature the description of a portfolio starts by defining the “losses” $`l_i`$ ($`i=1,\mathrm{},N`$) which takes values $`0`$ and $`1`$ depending on whether the borrower $`i`$ returns the credit or fails to return it. For convenience, we define $`S_i=12l_i`$ which takes values $`+1`$ when a borrower returns the credit and $`1`$ when it fails to return it. Note that the number of losses is $`L=_il_i=(N_iS_i)/2`$.

The simplest model is the one of independent borrowers. This is characterized by a single parameter which is the probability of default $`p`$. In this case, the loss probability is given by the binomial distribution:

$$p_B(L)=\left(\begin{array}{c}N\\ L\end{array}\right)(1p)^{NL}p^L.$$
(1)
The expected value of the number of losses is $`L=Nl_i=Np`$ and its variance is $`V(L)=N^2V(l_i)=N^2p(1p)`$. (The corresponding expressions in terms of the $`S_i`$ variables can be easily derived taking into account that $`S_i=N(12p)`$ and $`V(S_i)=4p(1p)`$.) The main drawback of this model is that it does not contain the possibility for correlations between borrowers which are known to exist and play an important role. Although in some cases the links between borrowers can be difficult to understand, in many cases borrowers tend to default in a cooperative manner.

The definition of the default correlation $`\rho `$ in terms of the $`S_i`$ variables is given by:

$$\rho =\frac{S_iS_jS_iS_j}{\sqrt{S_i^2S_i^2}\sqrt{S_j^2S_j^2}}(ij).$$
(2)
The same formula applies for the $`l_i`$ variables, as can be checked straightforwardly. To include correlations different theoretical models have been proposed by modifying the model of independent borrowers including phenomenological hypothesis on how $`p`$ depends on external factors R.Martin et al. (2001). For instance, the economic conditions can be described by a latent variable $`y`$ which is a random variable with a certain probability density $`f(y)`$. Then, the concept of default probability $`p`$ is substituted by the concept of conditional default probability $`\varphi (y)`$ such that for a fixed value of $`y`$ the probability of $`L`$ defaults is binomial with $`p=\varphi (y)`$. Therefore,

$$p(L)=_{\mathrm{}}^{\mathrm{}}𝑑yf(y)\left(\begin{array}{c}N\\ L\end{array}\right)\left[1\varphi (y)\right]^{NL}\varphi (y)^L.$$
(3)
Such conditionally independent models (CIM) are typically characterized by two parameters which could be empirically determined. First, the average (also called unconditional) default probability given by:

$$\overline{p}l_i=_{\mathrm{}}^{\mathrm{}}𝑑y\varphi (y)f(y).$$
(4)
Second, the default correlation which, for the case of homogeneous portfolios, is given by:

$$\rho =\frac{_{\mathrm{}}^{\mathrm{}}\varphi (y)^2f(y)𝑑y\left[_{\mathrm{}}^{\mathrm{}}\varphi (y)f(y)𝑑y\right]^2}{_{\mathrm{}}^{\mathrm{}}\varphi (y)f(y)𝑑y\left[_{\mathrm{}}^{\mathrm{}}\varphi (y)f(y)𝑑y\right]^2}.$$
(5)
For instance, the Merton based model (used by JPMorgan and KMV companies in their credit portfolio models) corresponds to the assumption that $`f(y)`$ is the $`N(0,1)`$ Gaussian probability density and the function $`\varphi `$ is:

$$\varphi _M(y)=\mathrm{\Phi }\left(\frac{c\sqrt{q}y}{\sqrt{1q}}\right)_{\mathrm{}}^{\frac{c\sqrt{q}y}{\sqrt{1q}}}\frac{1}{\sqrt{2\pi }}e^{z^2/2}𝑑z,$$
(6)
where $`\mathrm{\Phi }`$ is the error function, $`c`$ is the “critical value” and $`q`$ is the “asset correlation” (which should not be confused with the default correlation). For this case it is easy to check that the average default probability is $`\overline{p}=\mathrm{\Phi }(c)`$, and that the default correlation is:

$$\rho =_{\mathrm{}}^c_{\mathrm{}}^c\frac{dudv}{2\pi \sqrt{1q^2}}e^{(u^2+v^22quv)/2(1q^2)}.$$
(7)
Other examples of CIM are the CreditPortfolioView model (of McKinsey) and the CreditRisk$`+`$ model (of Credit Suisse Financial Products) which use different functions $`\varphi (y)`$ and $`f(y)`$ M.Crouhy et al. (2000); M.Gordy (2000). The CIM have been shown to only capture cyclic correlations, leaving aside direct correlations, contagion and cascade effects K.Giesecke and S.Weber (2002).

In this letter we propose a different method to model a credit portfolio that includes default correlations in a fundamentally different way. We describe portfolios with vectors $`X=(S_1,S_2,\mathrm{},S_N)`$. The set of all possible $`2^N`$ vectors $`X`$ is, therefore, our phase space $`\mathrm{\Omega }=\{X\}`$. To solve a credit portfolio model is to assign a probability law $`P(X)`$ on $`\mathrm{\Omega }`$. Instead of trying to derive which are the dependencies between the different borrowers or between the borrowers and the external factors, we will assume that default correlation $`\rho `$ is different from zero and derive $`P(X)`$ from the Maximum Statistical Entropy principle. This principle, which is essential in Information Theory, has been widely used for building probability distributions when a underlying theory is lacking. It states that one should consider the model $`P(X)`$ that maximizes the entropy functional:

$$S[P]=\underset{X\mathrm{\Omega }}{}P(X)\mathrm{ln}P(X),$$
(8)
subject to the conditions imposed by the previous known information. Such a method allows to control what are the exact ingredients (or the exact information) that is taken into account to formulate the model. This is opposite to the above presented theoretical models which assume a rigid structure for relations between borrowers that cannot be empirically checked.

It is instructive to first present the derivation of the simplest model. We will impose only two conditions. First the normalization condition:

$$\underset{X\mathrm{\Omega }}{}P(X)=1.$$
(9)
Second we assume that we have a previous knowledge of the expected value of the number of losses $`L`$ or, equivalently, of the average default probability $`\overline{p}`$. In terms of the probability law $`P(X)`$ this condition is written as:

$$S_i=\underset{X\mathrm{\Omega }}{}P(X)M(X)=N(12\overline{p}),$$
(10)
where we have defined the function $`M(X)=_iS_i`$. We want to stress that the knowledge of the average default probability $`\overline{p}`$ allows to fix the expected value of $`M`$ but not its exact value. Maximize (8) with conditions (9) and (10) is equivalent to maximize the Lagrange functional:

$`\mathrm{\Psi }_0[P]`$ $`=`$ $`{\displaystyle \underset{X\mathrm{\Omega }}{}}P(X)\mathrm{ln}P(X)\lambda \left({\displaystyle \underset{X\mathrm{\Omega }}{}}P(X)1\right)`$ (11)
$``$ $`\alpha \left({\displaystyle \underset{X\mathrm{\Omega }}{}}M(X)P(X)N(12\overline{p})\right),`$
where $`\lambda `$ and $`\alpha `$ are Lagrange multipliers. The result of the maximization is $`P(X)=e^{\lambda 1+\alpha M(X)}`$. The parameter $`\lambda `$ can be easily deduced imposing (9). To obtain $`\alpha `$ we require $`M(X)=N(12\overline{p})`$. A straightforward algebra renders $`\alpha =\mathrm{tanh}^1(12\overline{p})`$. The probability $`P(X)`$ results in:

$$P(X)=(1\overline{p})^{\frac{N+M(X)}{2}}\overline{p}^{\frac{NM(X)}{2}}.$$
(12)
Note that this probability depends on $`X`$ only through the function $`M(X)`$. Since there are $`\left(\genfrac{}{}{0pt}{}{N}{L}\right)`$ different vectors that give the same loss $`L`$, the probability $`p(L)`$ turns out to be Eq. (1) substituting $`p`$ by $`\overline{p}`$.

The inclusion of the existence of correlations between borrowers into the model requires to introduce a third condition in the maximization of (8). In this case, the symmetric function of $`X`$ whose expected value is fixed by the default correlation $`\rho `$ is:

$$H(X)=\frac{2}{N1}\underset{i>j}{}S_iS_j.$$
(13)
Note that the number of terms in the sum is $`N(N1)/2`$ so that the pre-factor $`2/(N1)`$ ensures that $`H(X)N`$. Using (2), (10) and the fact that $`S_i^2=1`$ it is straigforward to check that the knowledge of $`\rho `$ allows to fix the average value of $`H(X)`$ according to:

$$H(X)=N\left[14\overline{p}(1\overline{p})(1\rho )\right].$$
(14)
Imposing this condition, together with conditions (9) and (10), leads to the following Lagrange functional:

$$\mathrm{\Psi }[P]=\mathrm{\Psi }_0[P]\beta \left[\underset{X\mathrm{\Omega }}{}H(X)P(X)H(X)\right].$$
(15)
Diferentiating and imposing (9) to eliminate $`\lambda `$ one obtains:

$$P(X)=\frac{e^{\alpha M(X)+\beta H(X)}}{_{X\mathrm{\Omega }}e^{\alpha M(X)+\beta H(X)}}.$$
(16)
This can be be rewriten as $`P(X)=e^{}/_{X\mathrm{\Omega }}e^{}`$ with:

$$(X)=\frac{2\beta }{N1}\underset{i<j}{}S_iS_j\alpha \underset{i=1}{\overset{N}{}}S_i,$$
(17)
which is the canonical probability distribution corresponding to the well known long-range Ising model (LRIM) of a finite system with $`N`$ spins and a convenient redefinition of the exchange constant and the external field M.M.Tsypin and H.W.J.Blöte (2000). Note that $`P(X)`$ in (16), again, is only a function of $`M(X)`$ since $`H(X)=\frac{1}{N1}\left(M(X)^2N\right)`$. Therefore the probability of $`L`$ losses will be obtained by multiplying expression (16) by the binomial number $`\left(\genfrac{}{}{0pt}{}{N}{L}\right)`$.

To determine the values of $`\alpha `$ and $`\beta `$ in (16) that allow to impose conditions (10) and (14) for given values of $`\overline{p}`$ and $`\rho `$ is a difficult task that has been extensively studied in the field of Statistical Mechanics. Note that the problem when $`N\mathrm{}`$ can be easily solved by a saddle point method. This method renders a solution that coincides with the solution obtained by the mean field hypothesis H.E.Stanley (1983). But, such an approximation is of little interest here because of two reasons: credit portfolios, although being large, are quite far from the thermodynamic limit and such a mean field approximation inmediately leads to $`\rho =0`$ (except for the unrealistic case of $`\overline{p}=0.5`$).

Fig. 1 shows an example of the probability $`p(L)`$ derived by the LRIM corresponding to $`\overline{p}=0.02`$ and $`\rho =0.01`$ ($`\alpha =\mathrm{8.85\hspace{0.33em}10}^3`$, $`\beta =1.01148`$). The behaviour is compared with the binomial model (without correlation) and the Merton based model. In this last case the parameters ($`c=2.054`$ and $`q=0.0718`$) have been calibrated in order to have the same values of $`\overline{p}`$ and $`\rho =0.01`$. Note that the LRIM is the only model able to capture the existence of a second peak in the large $`L`$ region (negative $`M`$ region) that appears due to the existence of correlations which favour a collective behaviour of the borrowers.

Although the peak in the right(large $`L`$ region) is usually very small compared to the peak in the left (small $`L`$ region) it has strong implications in the evaluation of credit risk. There are different strategies to measure risk. A common one is the computation of the so called Value-at-Risk (VaR): this is the number of failures $`\mathrm{\Lambda }`$ that must be considered in order to ensure that the probability of having more than $`\mathrm{\Lambda }`$ fails is smaller than a certain tiny confidence level $`\alpha _c`$. Thus $`\mathrm{\Lambda }`$ is given by the solution of the equation:

$$\underset{L=\mathrm{\Lambda }}{\overset{N}{}}p(L)=\alpha _c$$
(18)
The smaller the confidence level, the larger $`\mathrm{\Lambda }`$. Fig. 2 shows an example of the behaviour of the sum in the left of (18) as a function of a variable $`\mathrm{\Lambda }`$. The horizontal line indicates the confidence level. Note that given the two-peaked shape of $`p(L)`$, the sum in (18) exhibitis a constant plateau. Consequently a small change in $`\rho `$ may produce a drastic change in the position of the VaR (crossing point).

From the physical point of view, this change is very similar to a first-order phase transition in which the order parameter (VaR) exhibits a discontinuity. The main difference is that our system is small and therefore real discontinuities do not exhist but, instead, there are sharp changes. In order to make clear this sharp behaviour we show in Figure 3 the VaR as a function of $`\rho `$ for $`\overline{p}=0.01`$ and $`\alpha _c=0.0001`$. The Merton model displays an smooth increase of the VAR with increasing correlation. For small correlations the differences between the LRIM and the binomial one are negligible. However, there exist a value of the correlations for which, according to the LRIM, a dramatic increase in the VaR occurs, which is not reproduced by the Merton model.

It is of paramount importance to locate the region where such strong change in the VaR occurs for different values of the confidence level $`\alpha _c`$. This is shown in the inset of Fig. 3. The linear axes represent the space of empirical parameters $`\overline{p}`$ and $`\rho `$. The different lines, corresponding to different confidence levels, separate the region in which the behaviour of the LRIM is similar to the binomial model (below) from the region where collective failure may occur and thus predictions with the binomial model strongly uderstimate the VaR. Other measurements of risk (like the expected shortfall) may exhibit even stronger diferences.

There are a series of possible extensions of the LRIM model: first, the extension to the case of non-homogeneous portfolios is straightforward. The borrowers need to be separated in different tranches. For each tranche the formulation of the model will require the knowledge of the default probability, the autocorrelation within each tranche and the correlations among different tranches. Note that in the LRIM model will be easy to introduce negative correlations between certain tranches to represent big competitors in small markets. This is impossible to reproduce with the conditional independent models. A second interesting extension will be to change the discrete variables $`S_i`$ to continuous variables $`\varphi _i`$. The Ising model, then, be substituted by a $`\varphi ^4`$ model. This will enable to describe the fact that when a borrower defaults there is always a recovery of part of the loan. A third extension could be the inclusion of randomness in the LRIM. For instance random fields and/or random bonds would allow to take into account the uncertainties in the knowledge of the empirical values of $`\rho `$ and $`\overline{p}`$. Finally one could include higher order interaction terms, like three spin terms. This would allow to introduce, in a controled way, the existence of three borrower correlations which are know to exist when loans are fully guaranteed by third parties. Note that the conditional independent models do also introduce higher order correlations, but in a completely uncontrolled manner.

Within the framework of information theory, it has been shown that the LRIM is the natural model to be used for the description of credit portfolios and risk modelling when knowledge about the existence of correlations is taken into account. One of the main results is that risk measurements can display big and sharp increases due to default correlations, in regions of the space parameters where other models fail to predict them. The discussion presented also shows that the LRIM can be useful not only in finance but, in general, for any complex system described by binary variables for which the available information is the average value and the correlations.

J.M. acknowledges fruitful discussions with Josep Masdemont. E.V. acknowledges fruitful discussions with A.Planes, J.Vives and R.S.Johal. This work has received financial support from CICyT (Spain), project MAT2001-3251 and CIRIT (Catalonia), project 2000SGR00025.


text:
# Tensor Product Variational Formulation for Quantum Systems

## 1 Introduction

As a precise numerical method for one-dimensional (1D) quantum systems, the density matrix renormalization group (DMRG) has been widely applied to unsolved problems in condensed matter physics.  If the linear dimension of the system is not so large ($``$ of the order of ten), the method is also applicable to 2D quantum models by use of a mapping from a finite size 2D cluster to a 1D chain that contains long-range interactions.

On the other hand, it is difficult to apply DMRG to infinitely large 2D quantum systems. Apparently the mapping from 2D lattice to 1D chain is inapplicable if the system size is infinite. In addition, the decay of the density matrix eigenvalue is very slow in higher dimensions,  and the phenomenon prevents to obtain a good renormalization group transformation that maps half (or quarter) infinite 2D system into a block spin. )))Existence of a good renormalization group transformation is not excluded, in case that the transformation is created some method other than the diagonalization of the density matrix. A naive extension of DMRG formulation to higher dimension shall encounter very poor numerical results, as was reported in an application of DMRG to 3D classical systems.

Stepping back to the DMRG formulation for 1D systems, the numerical efficiency of DMRG partially comes from its variational structure, where a trial state is represented as a product of orthogonal matrices.  It is possible to extend such a construction of variational state to 2D quantum and 3D classical systems. For these higher dimensional systems one has to prepare a trial states that are represented as 2D product of local weights. For example, Martín-Delgad, et al. employed the 6-vertex model as a trial state for 2D lattice spin/electron systems.  Okunishi and Nishino considered an extension of the Kramers-Wannier approximation  to the 3D Ising model, representing the trial state as the 2D Ising model under magnetic field,  or more generally as the interaction round a face (IRF) model.  Let us call these variational approaches as ‘the tensor product variational approach’ (TPVA) in the following.

In this paper we investigate the numerical efficiency of TPVA when it is applied to the square lattice $`S=1/2`$ Heisenberg and XY models. We employ an isotropic and uniform IRF model as the trial state, that contains 3 adjustable parameters. Our trial state includes the one parameter variational state proposed by Suzuki and Miyashita for the study of square lattice XY model.  The local construction of the variational state enables us to numerically calculate the energy expectation value by use of DMRG  or by Monte Carlo simulations.

In the next section we explain the construction of the trial state, and the way of calculating the energy expectation value. We show the numerical result in §3, and discuss a way of improving the variational formulation in the last section.

## 2 Construction of the Variational State

We consider the $`S=1/2`$ XXZ model on the square lattice as an example of 2D quantum systems. Its Hamiltonian is represented as

$$H=\underset{\{\mathrm{𝐫𝐫}^{}\}}{\overset{}{}}\left(S_𝐫^\mathrm{x}S_𝐫^{}^\mathrm{x}+S_𝐫^\mathrm{y}S_𝐫^{}^\mathrm{y}+\alpha S_𝐫^\mathrm{z}S_𝐫^{}^\mathrm{z}\right),$$
(2.1)
where $`\{\mathrm{𝐫𝐫}^{}\}`$ denotes pairs of neighboring lattice sites. The parameter $`\alpha `$ is chosen to be either 1 (the Heisenberg model) or zero (the XY model). Since the lattice is bipartite, the following Hamiltonian

$$\stackrel{~}{H}=\underset{\{\mathrm{𝐫𝐫}^{}\}}{\overset{}{}}\left(S_𝐫^\mathrm{x}S_𝐫^{}^\mathrm{x}S_𝐫^\mathrm{y}S_𝐫^{}^\mathrm{y}+\alpha S_𝐫^\mathrm{z}S_𝐫^{}^\mathrm{z}\right),$$
(2.2)
has the same same energy spectrum as $`H`$. In the following we treat $`\stackrel{~}{H}`$ instead of $`H`$ for the purpose of simplifying the variational formulation.)))One has to prepare two different local weights for $`H`$ and multiply them alternatively when constructing the variational state. Note that the mapping from $`H`$ to $`\stackrel{~}{H}`$ has nothing to do with the elimination of negative sign from the local weight. Our interest is in finding out a good variational function that minimizes the energy expectation value

$$\mathrm{\Lambda }=\frac{\mathrm{\Psi }|\stackrel{~}{H}|\mathrm{\Psi }}{\mathrm{\Psi }|\mathrm{\Psi }}$$
(2.3)
within the restriction where $`|\mathrm{\Psi }`$ is represented as a product of local weights.

Let us introduce a notation $`\sigma _𝐫^{}=2S_𝐫^z=\pm 1`$, and write the trial wave function as $`\mathrm{\Psi }(\{\sigma \})\{\sigma \}|\mathrm{\Psi }`$, where $`\{\sigma \}`$ represents a spin configuration of all the spins on the 2D lattice. We employ a trial wave function in the form of the uniform product of local weights

$$\mathrm{\Psi }(\{\sigma \})=\underset{𝐫}{\overset{}{}}W\left(\begin{array}{cc}\sigma _{𝐫+\widehat{𝐣}}& \sigma _{𝐫+\widehat{𝐢}+\widehat{𝐣}}\\ \sigma _𝐫^{}& \sigma _{𝐫+\widehat{𝐢}}\end{array}\right).$$
(2.4)
The weight $`W`$ is a function of 4 neighboring spin variables $`\sigma _𝐫^{}`$, $`\sigma _{𝐫+\widehat{𝐢}}`$, $`\sigma _{𝐫+\widehat{𝐣}}`$, and $`\sigma _{𝐫+\widehat{𝐢}+\widehat{𝐣}}`$ where $`\widehat{𝐢}`$ and $`\widehat{𝐣}`$ represent the unit lattice vector to X and Y directions, respectively. According to the symmetries of the Hamiltonian $`\stackrel{~}{H}`$, the local weight has only 3 independent parameters, that are

$`1`$ $`=`$ $`W\left(\begin{array}{cc}+& \\ & +\end{array}\right)=W\left(\begin{array}{cc}& +\\ +& \end{array}\right)`$ (2.9)
$`a`$ $`=`$ $`W\left(\begin{array}{cc}+& +\\ & \end{array}\right)=W\left(\begin{array}{cc}+& \\ +& \end{array}\right)=W\left(\begin{array}{cc}& \\ +& +\end{array}\right)=W\left(\begin{array}{cc}& +\\ & +\end{array}\right)`$ (2.18)
$`b`$ $`=`$ $`W\left(\begin{array}{cc}+& \\ & \end{array}\right)=W\left(\begin{array}{cc}& +\\ & \end{array}\right)=W\left(\begin{array}{cc}& \\ & +\end{array}\right)=W\left(\begin{array}{cc}& \\ +& \end{array}\right)`$ (2.27)
$`=`$ $`W\left(\begin{array}{cc}& +\\ +& +\end{array}\right)=W\left(\begin{array}{cc}+& \\ +& +\end{array}\right)=W\left(\begin{array}{cc}+& +\\ +& \end{array}\right)=W\left(\begin{array}{cc}+& +\\ & +\end{array}\right)`$ (2.36)
$`c`$ $`=`$ $`W\left(\begin{array}{cc}& \\ & \end{array}\right)=W\left(\begin{array}{cc}+& +\\ +& +\end{array}\right),`$ (2.41)
where we have written up and down spins by ’$`+`$’ and ’$``$’, respectively. Since we have constructed the trial state as a uniform product, the minimization of the variational ratio (Eq.(2$``$3)) is equivalent to the minimization of the local energy for an arbitrary bond

$$\lambda (a,b,c)=\frac{\mathrm{\Psi }|\left(S_𝐫^\mathrm{x}S_𝐫^{}^\mathrm{x}S_𝐫^\mathrm{y}S_𝐫^{}^\mathrm{y}+\alpha S_𝐫^\mathrm{z}S_𝐫^{}^\mathrm{z}\right)|\mathrm{\Psi }}{\mathrm{\Psi }|\mathrm{\Psi }},$$
(2.42)
where $`𝐫`$ and $`𝐫^{}`$ are neighboring lattice points. In addition, the product structure of the trial state enables us to obtain the denominator

$$\mathrm{\Psi }|\mathrm{\Psi }=\underset{𝐫}{\overset{}{}}\left[W\left(\begin{array}{cc}\sigma _{𝐫+\widehat{𝐣}}& \sigma _{𝐫+\widehat{𝐢}+\widehat{𝐣}}\\ \sigma _𝐫^{}& \sigma _{𝐫+\widehat{𝐢}}\end{array}\right)\right]_{}^2$$
(2.43)
as the partition function of the isotropic IRF model, which is specified by the local Boltzmann weight $`W_{}^2`$. Thus the norm of the trial state can be accurately calculated by use of DMRG applied to 2D classical systems.  Similarly the numerator is also a partition function of an IRF model that has additional structure around the bond $`\{\mathrm{𝐫𝐫}^{}\}`$ in Eq.(2$``$6), and can be calculated with sufficient numerical precision by DMRG as proposed by Hieida et al.  We use CTMRG,  which is a variant of DMRG, for the calculations of the following result.

## 3 Calculated Result

The minimum of the variational energy $`\lambda (a,b,c)`$ in Eq.(2$``$6) can be detected by way of the parameter sweep for $`a`$, $`b`$, and $`c`$. Table I shows the optimal parameter sets that give the lowest variational energy $`\lambda _{\mathrm{min}}^{}`$ for the isotropic Heisenberg model ($`\alpha =1`$) and the XY model ($`\alpha =0`$). The diagonal spin correlation $`S_𝐫^\mathrm{z}S_𝐫^{}^\mathrm{z}`$ and the off-diagonal one $`S_𝐫^\mathrm{x}S_𝐫^{}^\mathrm{x}`$ between the neighboring sites are also shown. In both cases the variational state is disordered, since there is no phase transition from the parameter limit $`a=b=c=1`$ to the shown parameter cases. Even at the Heisenberg point, the state is not isotropic, as is observed from the difference between $`S_𝐫^\mathrm{z}S_𝐫^{}^\mathrm{z}`$ and $`S_𝐫^\mathrm{x}S_𝐫^{}^\mathrm{x}`$.

Let us compare the calculated variational energy with two representative ground state energy estimations. Table II shows our result, a recent Monte Carlo (MC) result by Sandvik,  and the variational energy calculated by Suzuki and Miyashita in 1978.  Our result for the Heisenberg (or the XY) model is 2.3 % (or 1.2 %) higher than the MC result.

## 4 Discussions

As an application of DMRG for infinite size 2D quantum systems, we employ TPVA for the square lattice $`S=1/2`$ Heisenberg and XY models. The obtained variational state is disordered even for the Heiseberg point, the result which suggests the loss of the antiferromagnetic correlation. The problem may be improved by increasing the variational parameter, as Nishino et al have introduced auxiliary variables to the local weight when they applied TPVA to 3D classical models.

## Acknowledgments

The authors thank to K. Okunishi, Y. Hieida, Y. Akutsu, and S. Miyashita for valuable discussions. The present research is partially supported by a Grant-in-Aid from Ministry of Education, Science and Culture of Japan. A. G is supported by Japan Socienty for the Promotion of Science (P01192) and by VEGA grant No. 2/3118/23.