text:
# Computers with closed timelike curves can solve hard problems

## 1 Computing with closed timelike curves

The recent success in the field of quantum computation shows how the power of computation can be affected by the particular choice of physical model for a computer. By assuming a computer which operates according to the laws of quantum mechanics, Peter Shor was able to devise an algorithm to factor large numbers exponentially more efficiently than the best known classical algorithm .

This success leads one to ask: are there other physical models for computation which will also result in much more powerful algorithms? As long as one is speculating, one might as well speculate wildly; so let us consider computers with access to closed timelike curves (CTCs), which are thereby able to send information (such as the results of calculations) into their own past light cones .

I argue that such computers would be able to solve computationally difficult problems with amazing (indeed, almost magical) efficiency. In honor of Shor, I consider an algorithm for factoring large numbers; but it is easy to see that a very large class of computationally difficult problems, including NP-complete and PSPACE-complete problems, can be solved by the same trick.

For purposes of illustration, I assume that this computer has a particular register, timeRegister, which can be set from the future by messages sent via the CTC. The command used for setting timeRegister is timeSet(t,x), where t is the time at which to set the register and x is the value to which it should be set. In between resetting events, the register’s value remains unchanged.

Now consider the following program:

```
        input(N);
        timeRegister = -1;
        t = clock();
        p = timeRegister;
        if (p > 1) and (N mod p = 0) go to FINAL;
        p = 1;
        do
          p = p + 1;
        until (N mod p = 0) or (p > sqrt(N));
        if (p > sqrt(N)) p = N;
 FINAL  timeSet(t,p);
        output(p);
        end;

```

At the start of the program, the computer checks the value of timeRegister; if it divides N, the computer skips to the end of the program, sends the answer back in time, and outputs it. Otherwise the computer exhaustively searches until it finds a factor of N, then sends the factor back in time and outputs it.

What happens when this program is executed? If timeRegister does not divide N when it is checked, the computer searches until it finds a factor, then sends it back; so timeRegister will divide N when it is checked, which is a contradiction. On the other hand, if timeRegister does divide N when it is checked, it will skip to the end, send the factor back, and quit. While this situation is quite bizarre, it is self-consistent; and if only self-consistent evolutions can occur, then this must be what happens .

The operation of this computer is reminiscent of an old time-travel paradox. The brilliant young inventor receives a message from her future self, telling her that she is going to invent a time machine, and giving her the details of its construction. She duly builds the machine and demonstrates it. When she is old and famous, she sends a message back to her younger self, telling her that she is going to invent a time machine, and giving her the details of its construction.

This situation is self-consistent, but still very strange; the information on how to build a time machine appears out of nowhere. On the other hand, if the information hadn’t appeared, that would have been self-consistent as well; the young inventor need not have ever discovered the time machine.

In the computer program, by contrast, there is a contradiction if the information doesn’t appear. This situation is created by the inner search loop, which is guaranteed to find the answer sooner or later. Therefore this inner loop is necessary for the algorithm to function, even though the loop itself will never be executed.

## 2 A recursive version

There are a couple of loopholes in this argument that need to be addressed. The contradiction depends on the ability of the computer to execute the inner loop and find an answer. It is not very difficult, however, to choose N so large that it would take longer than the lifetime of the universe to find a factor by exhaustive search. It is quite easy to choose N big enough that it will take longer than the lifetime of any reasonable computer. If the program never gets to the final timeSet command, it will never send the factor back, and no contradiction can arise.

A similar problem might arise if the CTC doesn’t extend arbitrarily far into the future. A CTC with limited extent might only be able to send information back a short distance in time; if the inner loop only finishes executing in the future of the CTC, the program will be unable to send the factor back, and once again no contradiction arises.

This might seem an insuperable problem, but I believe that it can be circumvented by using a more clever algorithm than exhaustive search. Consider the following function:

```
        function factor(N,nStart,nEnd)
          timeRegister=-1;
          t = clock();
          p = timeRegister;
          if (p=0) or (N mod p = 0) go to FINAL;
          if (nEnd-nStart < nTractable)
            p = nStart-1;
            do
              p = p + 1;
            until (N mod p = 0) or (p > nEnd);
            if (p > nEnd) p = 0;
          else
            p = factor(nStart,(nStart+nEnd)/2);
            if (p = 0)
              p = factor(1+(nStart+nEnd)/2,nEnd);
 FINAL    timeSet(t,p);
          return p;
        end function;

```

The function factor(N,nStart,nEnd) looks for a factor of N within the range nStart to nEnd. If there is one, then it returns the factor; if not, it returns 0.

The function first checks to see if an answer has been sent back from the future. If one has, it skips to the end and returns it. If not, it checks if the range nEnd-nStart is small enough to search in a small time (determined by a constant parameter nTractable). If the range is small enough, it loops until it finds an answer. If not, it breaks the range into two parts and calls itself recursively for each subrange. At the end, it sends the answer back in time and returns it.

The factoring program then takes the following form:

```
        input(N);
        timeRegister = -1;
        t = clock();
        p = timeRegister;
        if (N mod p = 0) go to FINAL;
        p = factor(N,2,sqrt(N));
        if (p = 0) p = N;
 FINAL  timeSet(t,p);
        output(p);
        end;

```

Once again, at the beginning the program checks to see if the answer has been sent back from the future. If it has, then it skips to the end, sends the answer back in time, and outputs it. If not, then it enters the recursion. At the bottom level of the recursion is a loop that can be executed in a short time. The loop is only executed if the result of the loop is not sent back in time, but if the loop is executed then the result will be sent back in time. Therefore the loop will not be executed, and the answer will appear when checked. At the next higher level of recursion, the call to factor won’t be made, because the answer there will already have appeared; and so forth, all the way to the top of the recursion. The function factor will actually never be called at all. The only self-consistent outcome is that the program finds the correct answer when it checks timeRegister at the beginning.

The only requirement for this program to work is that the number of recursive calls to break the interval down into a tractable subinterval not be too big. Since the number of levels of recursion goes like $`\mathrm{log}_2𝙽`$, this is not very restrictive.

## 3 Harder problems

The particular algorithm I presented solved the factoring problem. While this problem is in NP, it is not NP-complete; but it is obvious that a program with the same structure could solve NP-complete problems as well. Indeed, it can solve even more difficult problems, as we shall see.

First, consider the satisfiability problem (SAT), which is to find a string of N bits $`x_1,\mathrm{},x_N`$ which simultaneously make true (satisfy) some set of clauses (i.e., logical statements) $`\varphi `$. These clauses $`\varphi `$ can be put together into a single logical statement involving the $`\{x_i\}`$, in conjunctive normal form. For example,

$$\varphi =(x_1\neg x_3x_{41})(\neg x_5x_{17})\mathrm{}.$$
(1)
The satisfiability problem can be solved (very inefficiently) by exhaustive search, merely trying every N-bit string until either finding one that satisfies the clauses, or determining that there isn’t one. By breaking down the set of all strings into smaller and smaller subsets using a recursive algorithm, one could modify the program in section II to solve SAT.

SAT is the canonical example of an NP-complete problem . Any problem in NP can be translated into an instance of the satisfiability problem with only polynomial overhead. Therefore, SAT is in that sense at least as difficult as any other problem in NP. However, these are not necessarily the most difficult problems that exist. Consider the following variant of the problem.

Once again we have a set of N bits and a set of clauses $`\varphi `$. This time, however, we don’t want to know if there is an assignment of bit values that satisfies $`\varphi `$; instead, we want to know if there is a value of $`x_1`$ such that for all $`x_2`$ there is a value of $`x_3`$ such that for all $`x_4`$, etc., such that $`phi`$ is satisfied:

$$x_1x_2x_3x_4\mathrm{}x_{N1}x_N\varphi ,$$
(2)
where I have assumed that N is even.

This is the quantified satisfiability problem (QSAT) . There is no known NP algorithm to solve QSAT; it is an example of a PSPACE-complete problem, i.e., a problem which can be solved using a computer with an amount of space polynomial in N, and which is polynomially equivalent to all other such problems. PSPACE-complete problems are believed to be strictly harder than NP-complete problems.

Such problems can be solved recursively. Suppose that $`\mathrm{𝚀𝚂𝙰𝚃}(\varphi )`$ is the function that evaluates QSAT for the set of clauses $`\varphi `$. Define $`\varphi _{00},\varphi _{01},\varphi _{10},\varphi _{11}`$, where $`\varphi _{ij}`$ is the set of clauses $`\varphi `$ with $`x_1x_2`$ replaced by $`ij`$. Then we see that

$$\mathrm{𝚀𝚂𝙰𝚃}(\varphi )=(\mathrm{𝚀𝚂𝙰𝚃}(\varphi _{00})\mathrm{𝚀𝚂𝙰𝚃}(\varphi _{01}))(\mathrm{𝚀𝚂𝙰𝚃}(\varphi _{10})\mathrm{𝚀𝚂𝙰𝚃}(\varphi _{11})).$$
(3)
The instances of QSAT on the right-hand side of (3) are all of length $`N2`$, and can be replaced by similar recursive expressions. By making use of this recursive structure, we can modify the program in section II to solve QSAT. Therefore, computers with CTCs should be able to solve not only NP-complete, but even PSPACE-complete problems efficiently.

One interesting difference in this case, however, is that while it is simple to check the answer to SAT (just by checking that the returned set of bit values does indeed satisfy $`\varphi `$), there is no efficient way of checking that an answer to QSAT is correct, in general.

## 4 Quantum considerations

So far in this paper I have treated both the computer and the CTC as if they were completely classical. I have adopted the assumption  that the allowed evolutions are those which do not produce a contradiction, and that the computer functions deterministically. How will these results change when we take quantum mechanics into account?

Of course, the most likely outcome of including quantum effects is that CTCs will no longer exist at all. This is the so-called “Chronology Protection Conjecture” of Steven Hawking, which postulates that the build-up of quantum fluctuations around a CTC will destabilize the spacetime and destroy the time machine . There is some evidence to support this conjecture, though there is (as yet) no proof.

Several authors have investigated quantum systems on a fixed background spacetime which includes CTCs. (See, e.g.,  and .) David Deutsch , in particular, has suggested that the Many-Worlds interpretation of quantum mechanics prevents time travel paradoxes. When one travels back in time to kill one’s grandfather (in the usual violent version of the paradox), one finds oneself in a different branch of the wavefunction; the future of the new “world” is changed, but not the old “world.”

Would this kind of argument eliminate a computer such as I describe? It is not obvious that it would. Deutsch’s argument prevents contradictions, but the operation of the algorithm, while mind-boggling, is not contradictory. The operation of the computer is deterministic, and should proceed identically in (almost) all universes in which it occurs (barring very improbable quantum fluctuations which, for example, demolish the lab). There seems no reason that self-consistent worlds with causal loops cannot exist. They don’t defy logic, but only common sense. Indeed, Deutsch himself in  suggested that closed timelike curves might make possible computers which solve hard problems.

It is possible that some other quantum effect might prevent the algorithm from working, while still allowing the existence of CTCs. But at present, no such argument has occurred to me.

## 5 Conclusions

It is very odd for information to suddenly appear out of nowhere, but in a universe with closed timelike curves such events can be expected to occur. It has widely argued that if CTCs are possible, the laws of physics should require that only internally consistent evolutions can occur, and that generalized versions of the principle of least action will enforce this behavior .

I’ve argued in this paper that one could exploit this tendency to design computers able to solve hard problems in very little time. I gave the specific example of factoring; but in section III, I argued that a similar algorithm could solve both NP-complete and PSPACE-complete problems as well, using the satisfiability and quantified satisfiability problems as examples. In all these cases, the answers appear out of nowhere in order to prevent logical contradictions from arising. Thus, these algorithms can be said to work because of the presence of brute-force search loops which are never actually executed.

This is a strange, though logically consistent, conclusion. But perhaps the best conclusion to draw is that it makes the existence of closed timelike curves even more unlikely.

## Acknowledgments

I would like to thank Steve Adler, Avi Wigderson, and Andrew Yao for reading drafts of this manuscripts and giving very helpful feedback; and Cara King, for encouraging me to write up this idea is spite of its innate craziness. This work was supported in part by the Martin A. and Helen Chooljian Membership in Natural Sciences, and by DOE Grant No. DE-FG02-90ER40542.


format: text

text:
# Transition from non Fermi Liquid Behavior to Landau Fermi Liquid Behavior Induced by Magnetic Fields

## Abstract

We show that a strongly correlated Fermi system with the fermion condensate, which exhibits strong deviations from Landau Fermi liquid behavior, is driven into the Landau Fermi liquid by applying a small magnetic field $`B`$ at temperature $`T=0`$. This field-induced Landau Fermi liquid behavior provides the constancy of the Kadowaki-Woods ratio. A re-entrance into the strongly correlated regime is observed if the magnetic field $`B`$ decreases to zero, then the effective mass $`M^{}`$ diverges as $`M^{}1/\sqrt{B}`$. At finite temperatures, the strongly correlated regime is restored at some temperature $`T^{}\sqrt{B}`$. This behavior is of general form and takes place in both three dimensional and two dimensional strongly correlated systems. We demonstrate that the observed $`1/\sqrt{B}`$ divergence of the effective mass and other specific features of heavy-fermion metals are accounted for by our consideration.

Recently, a peculiar critical point was observed in heavy-fermion metal YbRh2Si2 at low temperatures $`T`$ . This critical point is driven by the magnetic field $`B`$ which suppresses the antiferromagnetic order when reaches the critical value, $`B=B_{c0}`$, while the effective mass $`M^{}`$ diverges as $`M^{}1/\sqrt{BB_{c0}}`$ . The study of the magnetic field dependence of the coefficients $`A`$, $`\gamma _0`$, and $`\chi _0`$ in the resistivity, $`\rho (T)=\rho _0+\mathrm{\Delta }\rho `$, with $`\mathrm{\Delta }\rho =A(B)T^2`$, specific heat, $`C/T=\gamma _0(B)`$, and the magnetic $`ac`$ susceptibility, has revealed that YbRh2Si2 behaves as a true Landau Fermi liquid (LFL) for $`B>B_{c0}`$ and the well-known Kadowaki-Woods ratio $`A/\gamma _0`$  is preserved . In contrast, for $`B=0`$, YbRh2Si2 demonstrates a non Fermi liquid (NFL) behavior, and the resistivity follows a quasilinear temperature dependence $`\mathrm{\Delta }\rho T`$ down to 80 mK at which antiferromagnetic (AF) order takes place. At lower temperatures, the resistivity in AF-ordered state is described by $`\mathrm{\Delta }\rho T^2`$. A similar picture is observed in heavy-fermion compounds Ce$`M`$In5 ($`M`$=Ir, Co, and Rh), where the electronic specific heat $`C`$ shows more pronounced metallic behavior at sufficiently high magnetic fields . These observations are consistent with the de Haas van Alphen (dHvA) studies of CeIrIn5, which find that the effective mass decreases with increasing $`B`$ .

It is pertinent to note that heavy fermion metals are more likely three dimensional (3D) then two dimensional (2D). The origin of NFL behavior observed in heavy-fermion metals is still a subject of controversy . Moreover, the observed constancy of the Kadowaki-Woods ratio when $`BB_{c0}`$  leads to the failure of the standard model of heavy-fermion metals, when the mass renormalization is supposed to come from the exchange by soft magnetic fluctuations in a 2D spin fluid . As a result, we are left with even more complicated and challenging problems in the physics of strongly correlated electrons.

In this Letter, we study the nature of the critical behavior, assuming that the fermion condensation phase transition (FCQPT)  plays the role of the critical point. Analyzing the appearance of the fermion condensate (FC) which occurs beyond the point of FCQPT in an electron Fermi liquid and induces the NFL behavior, we show that the liquid is driven, by applying a weak magnetic field, back into a specific LFL state with effective mass $`M^{}1/\sqrt{B}`$. The LFL behavior induced by rather low magnetic fields, provides constancy of the Kadowaki-Woods ratio. But the strongly correlated regime is restored when the magnetic field $`B`$ decreases to zero, while the effective mass $`M^{}`$ diverges as $`M^{}1/\sqrt{B}`$. Also the strongly correlated regime is restored at some finite temperature $`T^{}(B)\sqrt{B}`$. Such a behavior is of general form and takes place in both three dimensional and two dimensional strongly correlated systems. We demonstrate that the observed crossover from NFL to LFL in certain heavy-fermion metals is accounted for by our consideration.

For the reader’s convenience we first outline the NFL behavior of Fermi systems with FC and the main properties of FCQPT using, as an example, a two-dimensional electron liquid in the superconducting state induced by FCQPT . At $`T=0`$, the ground state energy $`E_{gs}[\kappa (𝐩),n(𝐩)]`$ is a functional of the superconducting order parameter $`\kappa (𝐩)`$ and of the quasiparticle occupation function $`n(𝐩)`$ and is determined by the known equation of the weak-coupling theory of superconductivity (see e.g. )

$$E_{gs}=E[n(𝐩)]+\lambda _0V(𝐩_1,𝐩_2)\kappa (𝐩_1)\kappa ^{}(𝐩_2)\frac{d𝐩_1d𝐩_2}{(2\pi )^4}.$$
(1)
Here $`E[n(𝐩)]`$ is the ground-state energy of normal Fermi liquid, $`n(𝐩)=v^2(𝐩)`$ and $`\kappa (𝐩)=v(𝐩)\sqrt{1v^2(𝐩)}`$. It is assumed that the pairing interaction $`\lambda _0V(𝐩_1,𝐩_2)`$ is weak. Minimizing $`E_{gs}`$ with respect to $`\kappa (𝐩)`$ we obtain the equation connecting the single-particle energy $`\epsilon (𝐩)`$ to the superconducting gap $`\mathrm{\Delta }(𝐩)`$

$$\epsilon (𝐩)\mu =\mathrm{\Delta }(𝐩)\frac{12v^2(𝐩)}{2\kappa (𝐩)},$$
(2)
including the chemical potential $`\mu `$. Here the single-particle energy $`\epsilon (𝐩)`$ is determined by the Landau equation

$$\epsilon (𝐩)=\frac{\delta E[n(𝐩)]}{\delta n(𝐩)},$$
(3)
while the equation for the superconducting gap $`\mathrm{\Delta }(𝐩)`$ takes the form

$$\mathrm{\Delta }(𝐩)=\lambda _0V(𝐩,𝐩_1)\kappa (𝐩_1)\frac{d𝐩_1}{4\pi ^2}.$$
(4)
If $`\lambda _00`$, then, the maximum value of the superconducting gap $`\mathrm{\Delta }_10`$, and Eq. (2) reduces to that proposed in Ref.

$$\epsilon (𝐩)\mu =0,\mathrm{if}\kappa (p)0(0<n(p)<1)\mathrm{for}pL_{FC}:p_ipp_f.$$
(5)
At $`T=0`$, Eq. (5) defines a new state of Fermi liquid with FC, such that the modulus of the order parameter $`|\kappa (𝐩)|`$ has finite values in the FC range of momenta $`L_{FC}`$: $`p_ipp_f`$, while the superconducting gap can be infinitely small, $`\mathrm{\Delta }_10`$, in this range . Such a state can be considered as superconducting, with infinitely small value of $`\mathrm{\Delta }_1`$ so that the entropy of this state is zero. This state, created by the quantum phase transition, disappears at $`T>0`$. The FCQPT can be considered as a “pure” quantum phase transition because it cannot take place at finite temperatures. Generally, this quantum critical point should not represent the termination at $`T>0`$ of a line of continuous transitions. However it corresponds to a certain critical value of density $`x=x_{FC}`$ (the critical point of FCQPT) which is determined also by Eq. (5). Notice that at finite temperatures the FCQPT continues to have a strong impact on the system properties, up to a certain temperature $`T_f`$ above which FC effects become insignificant . FCQPT does not violate rotational or translational symmetry of the order parameter $`\kappa (𝐩)`$. It follows from Eq. (5) that the quasiparticle system “splits” into two quasiparticle subsystems: one in the $`L_{FC}`$ range, occupied by the quasiparticles with enhanced effective mass $`M_{FC}^{}1/\mathrm{\Delta }_1`$, and another with LFL effective mass $`M_L^{}`$ at $`p<p_i`$. If $`\lambda _00`$, then $`\mathrm{\Delta }_1`$ becomes finite, and the finite value of the effective mass $`M_{FC}^{}`$ in $`L_{FC}`$ can be obtained from Eq. (2) as

$$M_{FC}^{}p_F\frac{p_fp_i}{2\mathrm{\Delta }_1},$$
(6)
while the effective mass $`M_L^{}`$ is only weakly disturbed. Here $`p_F`$ is the Fermi momentum. It follows from Eq. (6) that the quasiparticle dispersion can be presented by two straight lines characterized by the effective masses $`M_{FC}^{}`$ and $`M_L^{}`$ respectively. These lines intersect near the electron binding energy $`E_0`$ which defines an intrinsic energy scale of the system:

$$E_0=\epsilon (𝐩_f)\epsilon (𝐩_i)\frac{(p_fp_i)p_F}{M_{FC}^{}}2\mathrm{\Delta }_1.$$
(7)
Let us assume that FC has just taken place, that is $`p_ip_Fp_f`$, the deviation $`\delta n(p)`$ from LFL occupation function is small (though finite), and $`\lambda _00`$. Expanding the functional $`E[n(p)]`$ in Eq. (3) in Taylor’s series with respect to $`\delta n(p)`$ and retaining the leading terms, we have

$$\mathrm{\Delta }E=\underset{\sigma _1}{}\epsilon _0(𝐩_\mathrm{𝟏},\sigma _1)\delta n(𝐩_\mathrm{𝟏},\sigma _1)\frac{d𝐩_1}{(2\pi )^2}+\underset{\sigma _1\sigma _2}{}F_L(𝐩_1,𝐩_2,\sigma _1,\sigma _2)\delta n(𝐩_\mathrm{𝟏},\sigma _1)\delta n(𝐩_\mathrm{𝟐},\sigma _2)\frac{d𝐩_1d𝐩_2}{(2\pi )^4},$$
(8)
where $`F_L(𝐩_1,𝐩_2,\sigma _1,\sigma _2)=\delta ^2E/\delta n(𝐩_1,\sigma _1)\delta n(𝐩_2,\sigma _2)`$ is the Landau interaction, and $`\sigma `$ denotes the spin states. Varying both sides of Eq. (8) with respect to the functions $`\delta n(p)`$ and taking into account Eq. (5), one obtains the FC equation

$$\mu =\epsilon (𝐩,\sigma )=\epsilon _0(𝐩,\sigma )+\underset{\sigma _1}{}F_L(𝐩,𝐩_1,\sigma ,\sigma _1)\delta n(𝐩_\mathrm{𝟏},\sigma _1)\frac{d𝐩_1}{(2\pi )^2};p_ipp_fL_{FC}.$$
(9)
Equation (9) acquires non-trivial solutions at the density $`x=x_{FC}`$ if the Landau amplitude $`F_L`$ (depending on density) is positive and sufficiently large, so that the potential energy integral on the right hand side of Eq. (9) prevails over the kinetic energy $`\epsilon _0(𝐩)`$ . Note, that in case of heavy fermion metals, this condition can be easily satisfied because of the huge effective mass. It is also seen from Eq. (9) that the FC quasiparticles form a collective state, since their energies are defined by the macroscopical number of quasiparticles within $`L_{FC}`$, and vice versa. The shape of their spectrum is not affected by the Landau interaction, which, generally speaking, depends on the system properties, including the collective states, impurities, etc. The only thing defined by the interaction is the width $`p_ip_f`$ of $`L_{FC}`$ (provided it exists). Thus, we can conclude that the spectra related to FC are of universal form.

At temperatures $`TT_c`$ when $`\mathrm{\Delta }_1`$ disappears, Eq. (6) for the effective mass $`M_{FC}^{}`$ is changed for

$$M_{FC}^{}p_F\frac{p_fp_i}{4T}.$$
(10)
The energy scale separating the slow dispersing low energy part, defined by the effective mass $`M_{FC}^{}`$, from the faster dispersing relatively high energy part, defined by the effective mass $`M_L^{}`$, can be estimated as $`E_0p_F(p_fp_i)/M_{FC}^{}`$ , so for the case of Eq. (10) it is

$$E_04T.$$
(11)
It follows from Eq. (10) that $`M_{FC}^{}`$ depends on the temperature, and the width $`\gamma `$ of the single-particle excitations results $`\gamma T`$, leading to a linear temperature dependence $`\mathrm{\Delta }\rho T`$ . This contrasts with the well-know LFL relations: $`\gamma T^2`$, and $`\mathrm{\Delta }\rho T^2`$.

Now we consider the behavior of an electronic system with FC in magnetic fields, supposing the coupling constant $`\lambda _00`$ to be infinitely small. As we have seen above, at $`T=0`$, the superconducting order parameter $`\kappa (𝐩)`$ is finite in the FC range, while the maximum value of the superconducting gap $`\mathrm{\Delta }_1\lambda _0`$ is infinitely small. Therefore, any small magnetic field $`B0`$ will destroy the coherence of $`\kappa (𝐩)`$ and thus FC itself. Also, the existence of FC can not be compatible with the evident Zeeman splitting of quasiparticle energy bands $`\epsilon (𝐩,\sigma )=\epsilon (𝐩)\sigma \mu _{eff}B`$ (see below for $`\mu _{eff}`$). To define the type of FC restructuring, simple energy reasons are invoked. On the one hand, the energy gain $`\mathrm{\Delta }E_B`$ due to the magnetic field $`B`$ is $`\mathrm{\Delta }E_BB^2`$ and tends to zero with $`B0`$. On the other hand, occupying the finite range $`L_{FC}`$ in the momentum space, FC delivers a finite gain in the ground state energy . Thus, a new state replacing FC should be very close in its ground state energy to the former state. Such a state is given by the multiconnected Fermi sphere, where the smooth quasiparticle distribution function $`n(𝐩)`$ in the $`L_{FC}`$ range is replaced by a multiconnected distribution $`\nu (𝐩)`$

$$\nu (𝐩)=\underset{k=1}{\overset{n}{}}\theta (pp_{2k1})\theta (p_{2k}p),$$
(12)
where the parameters $`p_ip_1<p_2<\mathrm{}<p_{2n}p_f`$ are adjusted to obey the normalization conditions:

$$_{p_{2k}}^{p_{2k+3}}\nu (𝐩)\frac{d𝐩}{(2\pi )^3}=_{p_{2k}}^{p_{2k+3}}n(𝐩)\frac{d𝐩}{(2\pi )^3}.$$
(13)
For the definiteness sake, we consider the most interesting case of 3D system, while the consideration of a 2D system goes along the same line. We note that the idea of multiconnected Fermi sphere, with production of new, interior segments of the Fermi surface, has been considered already . Let us assume that the thickness of each interior block is approximately the same $`p_{2k+1}p_{2k}\delta p`$ and $`\delta p`$ is defined by $`B`$. Then, the single-particle energy in the region $`L_{FC}`$ can be fitted by

$$\epsilon (𝐩)\mu \mu \frac{\delta p}{p_F}\left[\mathrm{sin}\left(\frac{p}{\delta p}\right)+b(p)\right].$$
(14)
The blocks are formed since all the single particle states around the minimum values of the fast sine function are occupied and those around its maximum values are empty, the average occupation being controlled by a slow function $`b(𝐩)\mathrm{cos}[\pi n(𝐩)]`$. It is seen from Eq. (14) that the effective mass $`m^{}`$ at each internal Fermi surface is of the order of the bare mass $`m_0`$, $`m^{}m_0`$. Upon substituting $`n(𝐩)`$ in Eq. (8) by $`\nu (𝐩)`$, defined by Eqs. (12) and (13), and taking into account the Simpson’s rule, we obtain that the minimum loss in the ground state energy due to the formation of blocks is about $`(\delta p)^4`$. This result can be understood if one bears in mind that the continuous FC function $`n(𝐩)`$ delivers the minimum value to the energy functional, Eq. (8), while the approximation $`\nu (𝐩)`$ by steps of size $`\delta p`$ produces the minimum error of the order of $`(\delta p)^4`$. On the other hand, this loss must be compensated by the energy gain due to the magnetic field. Thus, we arrive at

$$\delta p\sqrt{B}.$$
(15)
With an account taken of the Zeeman splitting in the dispersion law, Eq. (14), each of the blocks is polarized, since its outer areas are occupied only by spin-up quasiparticles. The width of this areas in the momentum space $`\delta p_0`$ is given by

$$\frac{p_F\delta p_0}{m^{}}B\mu _{eff},$$
(16)
where $`\mu _{eff}\mu _B`$ is the effective moment. We can consider such a polarization, without perturbing the previous estimates, since it is seen from Eq. (15) that $`\delta p_0/\delta p1`$. The total polarization $`\mathrm{\Delta }P`$ is obtained multiplying $`\delta p_0`$ by the number $`N`$ of the blocks which is proportional to $`1/\delta p`$, $`N(p_fp_i)/\delta p`$. Taking into account Eq. (15), we obtain

$$\mathrm{\Delta }Pm^{}\frac{p_fp_i}{\delta p}B\mu _{eff}\sqrt{B},$$
(17)
and thus it prevails over $`B`$ contribution from the LFL part. On the other hand, this quantity can be expressed as

$$\mathrm{\Delta }PM^{}B,$$
(18)
where $`M^{}`$ is the “average” effective mass related to the finite density of states at the Fermi level,

$$M^{}Nm^{}\frac{1}{\delta p}.$$
(19)
We can also conclude that $`M^{}`$ defines the specific heat.

Otherwise Eq. (15) can be examined, starting from a different point surmising that multiconnected Fermi sphere can be approximated by a single block. Let us put $`\lambda _0=0`$. Then, the energy gain due to the magnetic field $`\mathrm{\Delta }E_BB^2M^{}`$. The energy loss $`\mathrm{\Delta }E_{FC}`$ because of the restructuring of the FC state can be estimated by using the Landau formula which directly follows from Eqs. (8) and (9)

$$\mathrm{\Delta }E_{FC}=(\epsilon (𝐩)\mu )\delta n(𝐩)\frac{d𝐩^3}{(2\pi )^3}.$$
(20)
As we have seen above, the region occupied by variation $`\delta n(𝐩)`$ has the length $`\delta p`$, while $`(\epsilon (𝐩)\mu )(pp_F)p_F/M^{}`$. As a result, we have, $`\mathrm{\Delta }E_{FC}=\delta p^2/M^{}`$. Upon equating $`\mathrm{\Delta }E_B`$ and $`\mathrm{\Delta }E_{FC}`$ and taking into account Eq. (19), we arrive at the following equation

$$\frac{\delta p^2}{M^{}}\delta p^3\frac{B^2}{\delta p},$$
(21)
which coincides with Eq. (15).

It follows from Eqs. (17) and (18) that the effective mass $`M^{}`$ diverges as

$$M^{}\frac{1}{\sqrt{B}}.$$
(22)
Equation (22) shows that, by applying a magnetic field $`B`$, the system can be driven back into LFL with the effective mass $`M^{}(B)`$ dependent on the magnetic field. It was demonstrated that the constancy of the Kadowaki-Woods ratio is obeyed by systems in the strongly correlated regime when the effective mass is sufficiently large . Therefore, we are led to the conclusion that, by applying magnetic fields, the system is driven back into LFL where the constancy of the Kadowaki-Woods ratio is obeyed. Since the resistivity $`\mathrm{\Delta }\rho (M^{})^2`$ , we obtain from Eq. (22)

$$\mathrm{\Delta }\rho \frac{1}{B}.$$
(23)
At finite temperatures, the system persists to be LFL, but there is a temperature $`T^{}(B)`$ at which the polarized state is destroyed. To calculate the function $`T^{}(B)`$ , we observe that the effective mass $`M^{}`$ characterizing the single particle spectrum cannot be changed at $`T^{}(B)`$. In other words, at the crossover point, we have to compare the effective mass defined by $`T^{}`$, Eq. (10), and that defined by the magnetic field $`B`$

$$\frac{1}{M^{}}T^{}\sqrt{B}.$$
(24)
As a result, we obtain

$$T^{}(B)\sqrt{B}.$$
(25)
At temperatures $`TT^{}`$, the system comes back into the state with $`M^{}`$ defined by Eq. (10), and we observe the NFL behavior. It follows from Eq. (25), that a heavy fermion system at some temperature $`T`$ can be driven back into LFL by applying strong enough magnetic field $`BB_{cr}(T^{})^2`$. We can also conclude, that at finite temperature $`T`$, the effective mass of a heavy fermion system is relatively field independent at magnetic fields $`BB_{cr}`$ and show a more pronounced metallic behavior at $`BB_{cr}`$ since the effective mass is decreased, see Eq. (22). The same behavior of the effective mass can be observed in the dHvA measurements. We note that our consideration is valid up to temperatures $`TT_f`$.

Now we are in position to consider the nature of the field-induced quantum critical point in YbRh2Si2. The properties of this antiferromagnetic (AF) heavy fermion metal with the ordering Neèl temperature $`T_N=70`$ mK were recently investigated in Refs.. In AF state, this metal shows LFL behavior. As soon as the weak AF order is suppressed either by a tiny volume expansion or by temperature, pronounced deviations from LFL behavior are observed. The experimental facts show that the spin density wave picture is failed when dealing with the obtained data . We assume that the electron density in YbRh2Si2 is close to the critical value $`(x_{FC}x)/x_{FC}1`$ , so that this system can be easily driven across FCQPT. Then, in the AF state, the effective mass is given by Eq. (22) and the electron system of YbRh2Si2 possesses LFL behavior. When the AF state is suppressed at $`T>T_N`$ the system comes back into NFL. By tuning $`T_N0`$ at a critical field $`B=B_{c0}`$, the itinerant AF order is suppressed and replaced by spin fluctuations . Thus, we can expect absence of any long-ranged magnetic order in this state, and the situation corresponds to a paramagnetic system with strong correlation in the field $`B=0`$. As a result, the FC state is restored and we can observe NFL behavior at any temperatures in accordance with experimental facts . As soon as an excess magnetic field $`B>B_{c0}`$ is applied, the system is driven back into LFL. To describe the behavior of the effective mass, we can use Eq. (22) substituting $`B`$ by $`BB_{c0}`$

$$M^{}\frac{1}{\sqrt{BB_{c0}}}.$$
(26)
Equation (26) demonstrates the $`1/\sqrt{BB_{c0}}`$ divergence of the effective mass, and therefore the coefficients $`\gamma _0(B)`$ and $`\chi _0(B)`$ should have the same behavior. Meanwhile the coefficient $`A(B)`$ diverges as $`1/(BB_{c0})`$, being proportional to $`(M^{})^2`$ , and thus preserving the Kadowaki-Woods ratio, in agreement with experimental findings . To construct $`TB`$ phase diagram for YbRh2Si2 we use the same replacement $`BBB_{c0}`$ in Eq. (25)

$$T^{}(B)\sqrt{BB_{c0}}.$$
(27)
The phase diagram given by Eq. (27) is in good qualitative agreement with the experimental data . We recall that our consideration is valid at temperatures $`TT_f`$. The experimental phase diagram shows that the behavior $`T^{}\sqrt{BB_{c0}}`$ is observed up to $`150`$ mK  and allows us to estimate the magnitude of $`T_f`$ which can reach at least $`1`$ K in this system. It pertinent to note, that it follows directly from our consideration that the similar $`TB`$ phase diagram given by Eq. (27) can be observed at least in case of strongly overdoped high-temperature compounds. Except very close to the small values of both $`B`$ and $`T`$, because at $`TT_c`$ the magnetic field is to be $`B>B_c`$, here $`B_c`$ is the critical field suppressing the superconductivity.

To conclude, we have demonstrated that a new type of the quantum critical point observed in heavy-fermion metal YbRh2Si2 can be identified as FCQPT.

This work was supported in part by the Russian Foundation for Basic Research, No 01-02-17189.


text:
# Matrix models vs. Seiberg–Witten/Whitham theories

## 1.

It is well known that partition functions of matrix models are $`\tau `$-functions of integrable hierarchies of the Toda type . In the specific double scaling limit, these $`\tau `$-functions become $`\tau `$-functions of various reduction of the KP hierarchy . If one makes the simplest, large-$`N`$ (planar) limit, the partition function becomes the $`\tau `$-function of the dispersionless Toda hierarchy, which in turn becomes the $`\tau `$-function of the dispersionless (reductions of) KP hierarchy  after performing the continuum limit (which basically means working nearby a singularity of the partition function). All these dispersionless hierarchies are just Whitham equations over trivial solutions to integrable (Toda, KP) hierarchies.

When solving matrix models, most attention was paid to one-cut solutions where the limiting eigenvalue distribution spans one interval on the real axis . The results on multi-cut solutions  were few . Recently, Dijkgraaf and Vafa proposed  the new insight on the multi-cut large-$`N`$ limit of matrix models. Namely, they associated this limit with a Riemann surface and some related SW system. Its prepotential, which we prove here to be the logarithm of the large-$`N`$ partition function, is typically associated with the logarithm of some Whitham $`\tau `$-functions . This hints that the matrix matrix model in the large $`N`$ limit of multi-cut type describes the Whitham system over a non-trivial, finite-gap solution to integrable (Toda, KP) hierarchy. In particular, this solution passes to a finite-gap solution of (reductions of) the KP hierarchy in the continuum limit.

In this paper, we restrict ourselves with the simplest example of the Hermitean one-matrix model. We show that coefficients of the potential of the model gives rise to Whitham flows and manifestly construct this Whitham system. In fact, the authors of  associated the $`N=1`$ SUSY gauge theory studied in  with the SW system related to the multi-cut planar limit of matrix models. From the point of view of $`N=1`$ SUSY theory, these coefficients must be identified with couplings in the tree superpotential, while the SW moduli are associated with v.e.v.’s of the gluino condensates. This gives an interpretation of the results of  in the Whitham hierarchy terms.

## 2.

We call SW system 111Various properties of such systems can be found in . the following set of data222Our definition of the SW prepotential does not imply any connection with prepotentials of $`N=2`$ SUSY gauge theories . Moreover, the prepotentials discussed in this paper are rather related to superpotentials of $`N=1`$ SUSY theory .:

* a family $``$ of Riemann surfaces (complex curves) $`𝒞`$ whose dimesnion coincides with the genus;333This restriction can be waved, see examples in .
* a meromorphic differential $`dS`$ whose variations w.r.t. moduli of curves are holomorphic.

This data allows one to define the notion of prepotential  related to some integrable system .

Indeed, one can introduce variables

$$a_i_{A_i}𝑑S$$
(1)
where $`A_i`$ are $`A`$-cycles on $`𝒞`$. Then,

$$d\omega _i\frac{dS}{a_i}$$
(2)
are canonical holomorphic differentials on $`𝒞`$ (normalized so that $`_{A_i}𝑑\omega _j=\delta _{ij}`$). Then, introducing $`B`$-cycles conjugated to $`A`$-cycles: $`A_iB_j=\delta _{ij}`$, where $``$ means intersection, we obtain that

$$\frac{}{a_i}_{B_j}𝑑S=_{B_i}𝑑\omega _j=T_{ij}$$
(3)
is the period matrix of $`𝒞`$ and is therefore symmetric. Hence, there exists a prepotential $``$ such that

$$\frac{}{a_i}=_{B_i}𝑑S$$
(4)
## 3.

Let us consider the Hermitean one-matrix model. Its partition function is given by the integral over Hermitean $`N\times N`$ matrix $`M`$

$$Z_N=𝒟Me^{N\mathrm{tr}V(M)}$$
(5)
where $`𝒟M`$ is the Haar measure on Hermitean matrices and the potential $`V(x)`$ is a polynomial of degree $`n+1`$. After integrating out angular variables, we obtain

$$Z_N\underset{i}{}d\lambda _ie^{NV(\lambda _i)+_{ji}\mathrm{log}(\lambda _i\lambda _j)}$$
(6)
In the large $`N`$ limit, it is standard to introduce the density of eigenvalues

$$\rho (\lambda )\frac{1}{N}\underset{i}{}\delta (\lambda \lambda _i).$$
(7)
Then, (6) can be rewritten as

$$Z_N\underset{i}{}d\lambda _ie^{N^2\left[{\scriptscriptstyle \rho (\lambda )V(\lambda )𝑑\lambda }{\scriptscriptstyle \rho (\lambda )\rho (\lambda ^{})\mathrm{log}(\lambda \lambda ^{})𝑑\lambda 𝑑\lambda ^{}}\right]}$$
(8)
In the large $`N`$ limit, this integral can be evaluated by the saddle point method. We then assume $`\rho (\lambda )`$ to be a continuous function such that

$$\rho (\lambda )0,\rho (\lambda )𝑑\lambda =1$$
(9)
We then obtain the saddle point equation

$$V(\lambda )+\xi =2\rho (\lambda ^{})\mathrm{log}(\lambda \lambda ^{})𝑑\lambda ^{},\lambda ,\lambda ^{}\text{supp}(\rho )$$
(10)
where the support of the function $`\rho (\lambda )`$ comprises $`\lambda `$ such that $`\rho (\lambda )0`$. It emerges in this equation because of the first condition in (9). The constant $`\xi `$ in equation (10) is just the Lagrange multiplier for the second condition in (9).

However, in order to use analytic tools (the Cauchy problem), we must investigate not eq.(10) but its derivative

$$V^{}(\lambda )=2\overline{)}\frac{\rho (\lambda ^{})}{\lambda \lambda ^{}}𝑑\lambda ^{},\text{or}V^{}(\lambda )=\frac{\rho (\lambda ^{})}{\lambda \lambda ^{}}𝑑\lambda ^{}.$$
(11)
In order to solve this equation, we introduce the function

$$y(\lambda )\frac{2}{N}\underset{i}{}\frac{1}{\lambda _i\lambda }+V^{}(\lambda )=2\overline{)}\frac{\rho (\lambda ^{})}{\lambda ^{}\lambda }𝑑\lambda ^{}+V^{}(\lambda )$$
(12)
such that its imaginary part coincides with $`\rho (\lambda )`$ because of (11). In the large $`N`$ limit, it satisfies the equation

$$y^2(\lambda )V^2(\lambda )+4\frac{V^{}(\lambda )V^{}(\lambda ^{})}{\lambda \lambda ^{}}\rho (\lambda ^{})𝑑\lambda ^{}y^2(\lambda )V^2(\lambda )+f_{n1}(\lambda )=0$$
(13)
where $`f_{n1}(\lambda )`$ is a polynomial of degree $`n1`$. This means that the general solution to (11) is

$$y^2=V^2(\lambda )f_{n1}(\lambda )=\underset{i=1}{\overset{2n}{}}(\lambda \mu _i)$$
(14)
This equation describes a hyperelliptic curve of genus $`n1`$. It is, however, not arbitrary for a fixed potential $`V(\lambda )`$, because it follows from (13) that

$$y(\lambda )V^{}(\lambda )=W(\lambda )\frac{2}{\lambda }+O\left(\lambda ^2\right)|_\lambda \mathrm{},$$
(15)
i.e., $`\mu _i`$ are not independent. Here $`W(\lambda )`$ is the standard loop mean . One more restriction comes from the normalization condition in (9), and one is left with $`n1`$ free moduli. This is exactly what we need for SW system given on curve (14).

A solution to eq.(11) is parameterized by $`n1`$ moduli, i.e., these moduli span the moduli space of planar limits of the matrix model. The function $`\rho (\lambda )`$, which is imaginary part of $`y(\lambda )/(2\pi )`$, has the support on $`n`$ different branching cuts.

Note that within the standard matrix model framework, there are two more requirements that leave no moduli in solution. First of all, one can easily see that the sign of $`\rho (\lambda )`$ changes when coming to the next cut. This spoils non-negativity of $`\rho (\lambda )`$ and means that solution (14) is not stable. This means that the cuts with $`\rho (\lambda )`$ negative must shrink to produce double points. Thus, $`y(\lambda )`$ becomes proportional, besides the square root of a polynomial, to some other polynomial that has odd numbers of zeroes between cuts and, therefore, changes sign on every next cut.

We shall explain below that one can easily include these double points into the general SW and Whitham framework. Moreover, they allow one to construct more general Whitham systems.444V.Kazakov suggested to overcome non-stability of solutions without double points via some proper analityc continuation.

The second requirement looks more fundamental. Namely, returning to original equation (10), one has to check that the Lagrange multiplier $`\xi `$ is the same for every cut (while (11) only guarantees it is a constant on a cut). The difference of values of $`\xi `$ on two neighbour cuts is equal to

$$\xi _{i+1}\xi _i=_{\mu _{2i}}^{\mu _{2i+1}}y(\lambda )𝑑\lambda $$
(16)
where the integral runs from the right end of the left cut to the left end of the right cut. This gives $`n1`$ additional constraints and leaves no moduli (there still remains a freedom in the number of cuts). We must wave this requirement in order to make the Whitham system nontrivial. So, instead of just matrix models, it is better to speak about matrix-model-like solutions of Cauchy problem (11). So, at the moment we just ignore this last restriction and work with solutions to eq.(11), which we call matrix model solutions. We shall return to this point later.

## 4.

Now we associate an SW system with the planar limit of the matrix model. The family of genus $`g`$ curves is described by eq.(14) with the restriction (15). We describe it on the complex plane $`\lambda `$ by $`n=g+1`$ cuts, Fig.1. Besides canonically conjugated $`A`$\- and $`B`$-cycles, we also use the linear combination of $`B`$-cycles: $`\overline{B}_iB_iB_{i+1}`$, $`\overline{B}_{n1}B_{n1}`$. Therefore, $`\overline{B}`$-cycles encircle the nearest ends of two neighbour cuts, while all $`B`$-cycles goes from a given right end of the cut to the last, $`n`$-th cut. For the sake of definiteness, we order all points $`\mu _i`$ in accordance with their index so that $`\mu _i`$ is to the right of $`\mu _j`$ if $`i>j`$.

Fig. 1. Structure of cuts and contours.

The SW differential is

$$dS=y(\lambda )d\lambda $$
(17)
Its variations w.r.t. moduli is holomorphic on $`𝒞`$ (13) because all moduli are hidden in the polynomial $`f_{n1}(\lambda )`$:

$$\frac{dS}{\text{moduli}}=\frac{f_{n1}(\lambda )}{\text{moduli}}\frac{d\lambda }{y}$$
(18)
This expression is holomorhic, because the leading coeffient of $`f_{n1}(\lambda )`$ is fixed by the normalization condition (9), and the differentials $`\lambda ^k\frac{d\lambda }{y}`$ are holomorphic on the curve $`𝒞`$ (13) for $`k=0,1,\mathrm{},n2`$. Therefore, we introduce the variables

$$a_i=\frac{1}{2}_{A_i}y(\lambda )𝑑\lambda =\mathrm{Im}_{\mu _{2i1}}^{\mu _{2i}}y(\lambda )𝑑\lambda =_{\mu _{2i1}}^{\mu _{2i}}\rho (\lambda )𝑑\lambda ,i=1,\mathrm{},n1$$
(19)
that have meaning of “the occupation numbers” (or numbers of eigenvalues) associated with a given cut. Note also that

$$\frac{1}{2}_{A_n}y(\lambda )𝑑\lambda =\mathrm{Im}_{\mu _{2n1}}^{\mu _{2n}}y(\lambda )𝑑\lambda =1\underset{i=1}{\overset{n1}{}}a_i$$
(20)
which follows from (9). This is exactly the condition that fixes the leading coefficient of $`f_{n1}(\lambda )`$ and leaves $`n1`$ moduli. It means that

$$\frac{1}{2}_{A_n}𝑑\omega _i=\frac{1}{2}_{A_n}\frac{dS}{a_i}=1\text{ for all }i$$
(21)
Now one defines the prepotential

$$\frac{}{a_i}=_{B_i}𝑑S$$
(22)
This prepotenatial is equal to logarithm of the matrix model partition function, $`\mathrm{log}Z_N`$ in the planar limit. Indeed, let $`𝒟`$ be a set of contours $`A_1A_2\mathrm{}A_gA_n`$. The large $`N`$ partition function is555From now on, we consider symbols $``$ and $`res`$ with additional factors $`(2\pi i)^1`$ so that $`\mathrm{res}_0\frac{d\xi }{\xi }=\mathrm{res}_{\mathrm{}}\frac{d\xi }{\xi }=\frac{d\xi }{\xi }=1`$.

$$\mathrm{log}Z_N=\frac{1}{2}_𝒟y(\lambda )V(\lambda )𝑑\lambda +\frac{1}{4}_{𝒟\times 𝒟}y(\lambda )\mathrm{log}(\lambda \lambda ^{})y(\lambda ^{})𝑑\lambda 𝑑\lambda ^{}$$
(23)
We now calculate the derivative of $`\mathrm{log}Z_N`$ w.r.t. $`a_i`$:

$$\frac{\mathrm{log}Z}{a_i}=\frac{1}{2}_𝒟𝑑\lambda \frac{y(\lambda )}{a_i}\left(V(\lambda )_𝒟𝑑\lambda ^{}\mathrm{log}(\lambda \lambda ^{})y(\lambda ^{})\right)$$
(24)
The expression in the brackets on the rhs of (24) is a step function, which is equal to $`\xi _i`$ on each cut $`A_i`$ and its values on different cuts are (16)

$$\xi (\lambda )V(\lambda )_𝒟𝑑\lambda ^{}\mathrm{log}(\lambda \lambda ^{})y(\lambda ^{})=\{\begin{array}{c}\xi _1h_1\text{for}\lambda A_1,\hfill \\ \xi _1+_{\overline{B}_1}y(\lambda ^{})𝑑\lambda ^{}h_2\text{for}\lambda A_2,\hfill \\ \mathrm{}\hfill \\ \xi _1+_{\overline{B}_1\overline{B}_2\mathrm{}\overline{B}_{g1}}y(\lambda ^{})𝑑\lambda ^{}h_g\text{for}\lambda A_g,\hfill \\ \xi _1+_{\overline{B}_1\overline{B}_2\mathrm{}\overline{B}_{g1}\overline{B}_g}y(\lambda ^{})𝑑\lambda ^{}h_n\text{for}\lambda A_n.\hfill \end{array}$$
(25)
We therefore have

$$\frac{\mathrm{log}Z_N}{a_i}=_𝒟\frac{dS}{a_i}h(\lambda )=_𝒟\omega _ih(\lambda )=\xi _i+\xi _n=_{\overline{B}_i\overline{B}_{i+1}\mathrm{}\overline{B}_g}𝑑S_{B_i}𝑑S$$
(26)
and $`Z_N`$ in the planar limit can be, indeed, identified with the prepotential $`e^{}`$.

## 5.

One can learn two lessons from this fact. First of all, we can return to the interpretation of different $`\xi _i`$ on different cuts within matrix model. The standard matrix model case of equal $`\xi _i`$’s can be now formulated as the set of conditions

$$\frac{}{a_1}=\mathrm{}=\frac{}{a_g}=0$$
(27)
These are the conditions of minimum of the matrix model partition function w.r.t. the occupation numbers. They can be removed by introducing different chemical potentials for different cuts.666Putting differently, one can interpret these conditions as a criterium of stability against tunneling of eigenvalues between different cuts . Stability is achieved by imposing equality of the chemical potentials of all cuts. We do not enter here any further details and go instead to another lesson.

We know from studies of matrix models that their partition functions are $`\tau `$-functions of some integrable hierarchies . What are they in the planar limit? We have just proved that such a partition function is an SW prepotential in this limit. One typically associates logarithms of Whitham $`\tau `$-functions with SW prepotentials . Therefore, we may expect that the matrix model partition function becomes the $`\tau `$-function of some Whitham hierarchy. An additional evidence for this comes from looking at the simplest one-cut large-$`N`$ solution of the matrix model, when the partition function becomes the $`\tau `$-function of the dispersionless Whitham hierarchy . Now we construct this hierarchy in very manifest terms.

First, we return to the problem of double points. Let us assume that some of the cuts shrink, i.e.,

$$y(\lambda )=M_{nk}(\lambda )\sqrt{\underset{i=1}{\overset{2k}{}}\left(\lambda \mu _i\right)}M_{nk}(\lambda )\sqrt{g_{2k}(\lambda )}$$
(28)
where $`M_{nk}(\lambda )`$ is a polynomial of degree $`nk`$ and $`g_{2k}(\lambda )`$ is a polynomial of degree $`2k`$. This means that one is effectively left with a new curve

$$y(\lambda )=\sqrt{g_{2k}(\lambda )}$$
(29)
This curve of lower genus $`k1`$ along with the differential $`dS=M_{nk}(\lambda )y(\lambda )d\lambda `$ remarkably give rise to a new SW system that depends on $`k1`$ moduli.

To see this, one needs to take into account that there still holds eq.(13),

$$y^2(\lambda )M_{nk}^2(\lambda )=V^2(\lambda )f_{n1}(\lambda )V^2(\lambda )2(V^{}(\lambda )W(\lambda ))_+,$$
(30)
where we let $`()_+`$ denote the polynomial part of the expression in brackets. Then, varying $`dS`$ and using (29), we obtain for the general variation $`\delta dS`$:

$$\delta dS=\delta \left(M_{nk}(\lambda )y(\lambda )\right)d\lambda =\frac{g_{2k}(\lambda )\delta M_{nk}(\lambda )+\frac{1}{2}M_{nk}(\lambda )\delta g_{2k}(\lambda )}{y}d\lambda $$
(31)
On the other hand, doing a variation $`\stackrel{~}{\delta }`$ of $`M_{nk}(\lambda )y(\lambda )`$ that does not alter the potential, we obtain from (30) that

$$\stackrel{~}{\delta }dS=\frac{1}{2}\frac{\stackrel{~}{\delta }f_{n1}(\lambda )}{M_{nk}(\lambda )y(\lambda )}d\lambda .$$
(32)
Because this variation is a particular case of (31), we obtain that zeroes of $`M_{nk}(\lambda )`$ in the denominator of (32) must cancel, so the maximum degree of the polynomial in the numerator is $`n2`$. The variation is then holomorphic on curve (29).

This solves the problem of double points. The corresponding system with double points (the large-$`N`$ limit of the matrix model) is still described by the SW theory.

## 6.

Let us return to the case $`n=k`$. This SW system is described by $`n1`$ and $`n+1`$ additional parameters. They are expected to be Whitham times giving flows on the moduli of the finite gap solution. In this moduli space there are deformations leaving the curve within the family (variations of $`f_{n1}(\lambda )`$ that does not change the potential), and those transversal to the family. These latter are defined by the potential. If one stays within the family of hyperelliptic curves (13), there are exactly $`n+1`$ transversal deformations. Therefore, in order to have enough many deformations (Whitham times), one needs to involve potentials of high enough degree, i.e. to deal with the construction with double points.

In our manifest construction of the Whitham system we mainly follow  (see also ). In order to construct a Whitham system, one needs to add to the SW data a set of punctures with local coordinates in their vicinity. These points here are the two infinities on the curve (29) and the local parameter is $`\eta =\frac{1}{\lambda }`$. Now one introduces a set of meromorphic differentials $`d\mathrm{\Omega }_n`$ with the poles only at punctures (since the hyperelliptic curve (29) is invariant w.r.t. the involution $`yy`$, from now on we just work with either of the two infinities, see ) and the behaviour

$$d\mathrm{\Omega }_m=\left(\eta ^{m1}+O(1)\right)d\eta ,\eta 0$$
(33)
Then, the Whitham system is generated by a set of equations for these differentials and the holomorphic differentials $`d\omega _i`$:

$$\frac{d\mathrm{\Omega }_p}{t_m}=\frac{d\mathrm{\Omega }_m}{t_p},\frac{d\mathrm{\Omega }_m}{a_i}=\frac{d\omega _i}{t_m},\frac{d\omega _i}{a_j}=\frac{d\omega _j}{a_i}$$
(34)
These equations implies that there exists a differential $`dS`$ such that

$$\frac{dS}{a_i}=d\omega _i,\frac{dS}{t_m}=d\mathrm{\Omega }_m$$
(35)
Let us check that the differential $`dSM_{nk}(\lambda )y(\lambda )d\lambda `$ given on the curve (29) with the relation for moduli (30) really satisfies (35).

Indeed, we have proved the first set of relations (35) in the previous paragraph. Now let us consider variations of the potential, i.e., variations w.r.t. Whitham times $`t_m`$. Then, we obtain instead of (32)

$$\delta dS=\frac{1}{2}\frac{\delta \left(V^2(\lambda )f_{n1}(\lambda )\right)}{M_{nk}(\lambda )y(\lambda )}d\lambda $$
(36)
while (31) still holds. Repeating the argument of the previous paragraph, we conclude that the zeroes of $`M_{nk}(\lambda )`$ cancel from the denominator and, therefore, the variation may have pole only at $`\lambda =\mathrm{}`$ or $`\eta =0`$, i.e. at the puncture. In order to estimate this pole, one needs to use (31), which implies that $`dS=M_{nk}(\lambda )y(\lambda )d\lambda (V^{}(\lambda )+O(\frac{1}{\lambda }))d\lambda `$ and, therefore, the variation of $`dS`$ at large $`\lambda `$ is completely determined by the variation of $`V^{}(\lambda )`$. Parameterizing $`V(\lambda )=^{n+1}t_m\frac{\lambda ^m}{m}`$ one comes to (35) up to a linear combination of holomorphic differentials. One may fix the normalization of $`d\mathrm{\Omega }_m`$ that are also defined up to a linear combination of holomorphic differentials so that eq.(35) would be exact. What does this normalization mean? Throughout all our consideration we deal with $`a_i`$’s and $`t_k`$’s as independent variables. This unambiguously defines the way $`a_i`$’s depend on the coefficients of $`f_{n1}`$ and is achieved merely by imposing the (obvious) condition

$$\frac{a_i}{t_m}=_{A_i}𝑑\mathrm{\Omega }_m=0i,m,$$
(37)
Thus, similarly to eq.(1) we can invariantly introduce variables $`t_m`$ via the relation

$$t_m=\mathrm{res}_{\eta =0}\eta ^mdS$$
(38)
and define the prepotential that depends on both $`a_i`$ and $`t_m`$ via the old relation (4) and the similar relation

$$\frac{}{t_m}=\frac{1}{m}\mathrm{res}_{\eta =0}\eta ^mdS$$
(39)
One can immediately prove that such a prepotential exists , i.e., the second derivatives are symmetric, and, moreover, similarly to Sec. 4, we find that thus defined $``$ coincides with $`\mathrm{log}Z_N`$ in the planar limit. For this, we apply the formula similar to (24) with the only difference that the potential $`V(\lambda )`$ itself is changed. We then obtain from (25)

$$\frac{\mathrm{log}Z}{t_m}=\frac{1}{2}_𝒟𝑑\lambda \frac{y(\lambda )}{t_m}\xi (\lambda )\frac{1}{2}_𝒟𝑑\lambda y(\lambda )\frac{\lambda ^m}{m}=\underset{i=1}{\overset{g}{}}\frac{a_i}{t_m}(h_ih_n)+\frac{1}{m}\mathrm{res}_{\eta =0}\eta ^mdS,$$
(40)
which by virtue of (37) gives (39). Therefore, $`Z_N`$ in the planar limit is the Whitham $`\tau `$-function, and the whole machinery of Whitham systems works here in full strength.

## 7.

After having constructed the large $`N`$ (planar) limit, next step is to take the (double scaling) continuum limit. Namely, one has to dwell nearby a singularity (branching point) of $`\rho (\lambda )`$. Say, one can work nearby the left end of the very right, $`n`$-th cut $`\mu _{2n1}`$ . The standard argument then is that one feels no other cuts, since they are far away. This would mean that multi-cut solutions coincide with the one cut solution in the continuum limit. This is, however, the case only if all the other branching points do not come close to $`\mu _{2n1}`$. Otherwise, there exist non-trivial continuum limits . In the most non-trivial situation, all the branching points but $`\mu _{2n}`$ come close to each other777Note that the authors of the paper  considered the branching point $`\mu _{2n1}`$ also located far away, at some large distance $`\mathrm{\Lambda }`$. Then, they were interested in the pieces of the prepotential that shows up a logarithmic behaviour at $`\mathrm{\Lambda }\mathrm{}`$ plus terms constant in $`\mathrm{\Lambda }`$. This degenerated situation better suits the $`N=1`$ SUSY theory needs, where $`\mathrm{\Lambda }`$ plays a role of the $`\mathrm{\Lambda }_{QCD}`$ parameter . In order to get this limit, one just suffices to rescale $`\rho (\lambda )\frac{\rho (\lambda )}{\sqrt{\lambda \mathrm{\Lambda }}}`$ with all the $`B`$-cycles still encircling the point $`\mathrm{\Lambda }`$. Thus, the number of cuts becomes $`n1`$, but still there is a puncture at infinity at $`\mathrm{\Lambda }`$. With properly substracted $`\mathrm{\Lambda }`$-divergent pieces, it gives the same function $`(a)`$, with $`a_i`$ just rescaled. This gives the calculational recipee of the paper . The only subtlety is that, in order to properly cut-off logarithmic terms, more concretely, to reproduce the terms $`a_i^2\mathrm{log}\mathrm{\Lambda }`$ in the prepotential, one needs to rescale properly the matrix model partition function by volumes of the unitary groups, see . For instance, the simplest non-trivial example of the two-cut solution
$$y=\sqrt{(\lambda ^2\mu ^2)(\lambda \mathrm{\Lambda })}$$
(41) reduces in this limit to the semi-circle distribution $`y=\sqrt{\lambda ^2\mu ^2}`$, but the prepotential is determined by the integral $`\frac{}{a}=_\mu ^\mathrm{\Lambda }\sqrt{\lambda ^2\mu ^2}𝑑\lambda `$ with $`a=_\mu ^\mu \sqrt{\lambda ^2\mu ^2}𝑑\lambda `$. This immediately gives $`a^2\mathrm{log}\frac{a}{\mathrm{\Lambda }^2}+\mathrm{}`$ and coincides with the result obtained by exact calculating with (41) and then bringing $`\mathrm{\Lambda }`$ to infinity, although $`a`$ gets rescaled by $`\sqrt{\mathrm{\Lambda }}`$. . This is equivalent just to sending $`\lambda _{2n}`$ to infinity. Such a curve describes a finite gap solution to the KdV hierarchy, moreover, the corresponding SW system is also associated with KdV . Therefore, we expect that the matrix model partition function in this limit describes (in leading order) the Whitham hierarchy over KdV finite gap solution.

However, it would be very instructive to construct an entire double scaling limit in this situation, in particular, to fix proper scaling behaviours. This means to match properly the growth of $`N`$ and approaching to the singularity. Then, one could address the problem of exact (matrix model)$``$ (SW) correspondence, in particular, in integrable terms. In particular, it would be interesting to see what is the proper deformation of SW systems in this case.

In this respect, the problem of studying higher-genus corrections looks very natural because the whole matrix-model-like solution (the solution to the loop equation, see ) in all genera is completely determined by the set of data $`\{a_i,t_m\}`$ and because it was proved  that spectral correlators manifest the universality property for multi-cut solutions as well. The integrable system that will appear in this approach must be a generalization of a Whitham system.

At last, let us note that the construction considered in this paper is directly extendable to other matrix models. In particular, one can consider the model of normal matrix that has much to do with the problem of Laplacian growth . In fact, its naive large $`N`$ limit describes how external and internal moments of a domain are related. This domain is a counterpart of the one cut. Moreover, the system is described by the Whitham hierarchy that is the dispersionless Toda system. Now, considering a multi-domain solution and introducing chemical potentials for different domains, one has to get the Whitham system over a finite-gap solution to the Toda system. This Whitham hierarchy should relate the external and internal moments of several domains.

The authors are grateful to A.Gorsky, S.Gukov, V.Kazakov, A.Marshakov and especially to A.Morozov for fruitful discussions. We are also due to Chung-I Tan for providing us with proper references on the multi-cut solutions.

We acknowledge the remarkable atmosphere at Mologa meeting where this work has been started.


text:
# 1 Introduction

## 1 Introduction

Germinal centers (GC) are an important part of the humoral immune response . They generate new plasma cells that produce antibodies of high affinity to a specific antigen. This affinity maturation process allows the organism to adapt its antibody repertoire to new antigens, to form the memory B cell compartment, and to react to antigens more efficiently compared to an immune response using B-cells from the repertoire only.

Typically, GCs are found in lymph nodes or in the spleen. In reaction to primary immunization some B-cells from the repertoire are activated and become sensible to proliferation progression factors. These B blasts (centroblasts) migrate towards the primary follicle system. In the environment of follicular dendritic cells and T cells they start a process of monoclonal expansion . After about three days the centroblast population reaches about 12000 cells due to an extremely short cell cycle time of $`6`$-$`7`$ hours . During this early phase of the GC reaction a reduction of diversity has been observed where those B cells of higher affinity to the antigen have been preserved . This has been interpreted as an early antigen-dependent preselection process of un-mutated B cells. The phase of monoclonal expansion is believed to be followed by a phase of frequent somatic hypermutation , giving rise to a great variety of different antibodies that are encoded in the blasts. These cells start to differentiate to centrocytes that are no longer in cell cycle but express large numbers of antibodies . Centrocytes have initiated a process of apoptosis . They are rescued from apoptosis in dependence on an interaction with antigens that are presented on the FDCs . All centrocytes that are not selected by such an interaction process die – and this is the large majority. Only positively selected centrocytes get the chance to interact with T cells and, if again selected, to further differentiate into antibody producing plasma cells or memory cells . As this selection process is believed to depend on the affinity of antibodies and antigen, the GC reaction provides an antibody optimization procedure.

GCs show some very interesting properties. Among those are the separation of centroblasts and centrocytes into two zones, the dark and the light zone . The origin of those zones has been discussed before in the framework of a stochastic space-time model for the GC reaction . Another property is the nearly all-or-none behavior of GCs . GC reactions of moderate intensity have rarely been observed. The reaction is either successful (concerning affinity maturation and number of produced output cells) or strongly suppressed. The reason for this behavior has not been resolved until now. Especially, the all-or-none behavior could not be reproduced in a previous modeling work .

One suggestion goes back to the starting point of the GC reaction, the originally activated B cells . The authors speculate that the mouse colony in their (pathogen-free) experiment may undergo a qualitatively and quantitatively diverse background immune response, i.e. leading to a diversification of the B cell repertoire between individuals. The consequence would be that the seeder cells of the GC reaction would have different affinities to the antigen, implying different GC reactions. This leads us to the hypothesis that the quality of GC seeder cells determines the fate of the GC reaction  and, in addition, may imply the all-or-none behavior of GCs. In the present article this hypothesis is discussed in the framework of a previously introduced stochastic space-time model for the GC reaction .

## 2 Methods

The model providing the basis of the analysis has been introduced in all details before . Therefore, only those model elements are pointed out in details that are of direct relevance in the present context. The basic concepts are only shortly summarized in the following.

Cells are represented on a fixed, equidistant, cubic, and two dimensional lattice. Each cell occupies exactly one node of the lattice. All cells can move on the lattice with some velocity, which is anti-proportional to the cell radius. The movement is undirectioned and restricted by the space available on the lattice. Each cell is of specific type (centroblast or centrocyte) which determines their possible actions: Centroblasts can divide. Each division allows a somatic hypermutation to occur with probability $`0.5`$. They also can differentiate to centrocytes provided that the differentiation process has been activated before. This happens if the centroblast binds a certain amount of signal molecules, which is secerned by the FDCs and diffuses over the lattice. Such a non-local interaction process turned out to be crucial for the establishment of the dark zone . Centrocytes die by apoptosis and are then eliminated from the lattice. They may be rescued by an interaction with antigens on FDCs, which becomes possible when centrocytes and FDCs are in direct contact on the lattice and provided that the affinity of the presented antibodies to the antigen is strong enough. Rescued centrocytes can further differentiate to plasma- or memory-cells (output cells) or restart to proliferate (i.e. become centroblasts again). Output cells leave the GC by diffusion. The parameters that have been used here and that have been determined before  are summarized in Tab. 1.

Most important is the description of affinity between antigen and antibodies. The model uses the traditional shape space concept . The dimension of the shape space is thought to represent phenotype properties of the antibodies. The affinity to the antigen reaches a maximum at the shape space position corresponding to optimal complementarity. A mutation (only those of phenotypical relevance are counted) is represented by a jump to a neighbor point in the shape space. The principle of affinity neighborhood  ensures that the affinity may be increased by a sequence of somatic hypermutations. One may think of climbing up an affinity hill. However, key mutations, that has been frequently discussed (especially for immunization with (4-hydroxy-3-nitrophenyl)acetyl (NP) ) are not easily covered by this concept. Therefore, the outcome of the analysis will have to be interpreted in this context.

The quality of a B cell is defined by its distance to the antigen (i.e. to the optimal B cell clone). The distance is calculated by the number of (phenotypically relevant) somatic hypermutations that are necessary to reach the antigen position in the shape space. The number of mutations is denoted by mutation-distance in the following. This corresponds to the definition of a 1-norm on the shape space which is assumed to have four dimensions. Using this norm the affinity of antibody and antigen is quantitatively described by a Gaussian function which is centered at the antigen position $`\mathrm{\Phi }_0`$ in the shape space . In the following this concept is used to describe the efficiency of the GC reaction in a qualitative and quantitative way. At first, the quality of a GC reaction is measured by the success of the affinity maturation process. To this end the fraction of high affinity output cells (the sum of plasma and memory cells) resulting from the GC reaction is analyzed. High affinity means that the antibody binds the antigen with a probability of more than $`30\%`$. Secondly, the number of resulting output cells (summed over the whole GC reaction) is a good measure for the quantitative efficiency of the GC reaction.

However, the above notion of B cell quality hides an intrinsic ambiguity. Let the antigen be at position $`\mathrm{\Phi }_0=(3,3,3,3)`$ and let the seeder cells have a mutation-distance of $`5`$. The outcome of the GC reaction depends on the exact position of the seeder cell, which may have the initial position $`\mathrm{\Phi }_1=(8,3,3,3)`$ or $`\mathrm{\Phi }_2=(6,5,4,4)`$. Assuming the mutations to occur randomly , the probabilistic weights $`p(\mathrm{\Phi })`$ of both seeder cells to find the antigen position $`\mathrm{\Phi }_0`$ considerably differs: $`p(\mathrm{\Phi }_1)=1/256`$ and $`p(\mathrm{\Phi }_2)=90/256`$. The difference is due to the number of possible permutations in the sequence of mutations. In other words, there exist more combinatorial possibilities for the mutations that all lead to the same target $`\mathrm{\Phi }_0`$. This ambiguity is mirrored in the GC efficiency (quality and quantity of output cells, see Fig. 1). Less probable mutation paths induce GC reactions of smaller efficiency. However, at least the fraction of high affinity output cells saturates for more probable mutation paths. In order to ensure a well-defined mutation distance, only seeder cells that reach the antigen on a mutation path of maximum probability will be taken into account. The results have been cross-checked using mutation paths of moderate probability.

Using the thus defined model concepts we initiate GC reactions with various types of seeder cells. The seeder cells are placed at random positions in the FDC network. The mutation-distance is varied between $`0`$ and $`12`$. Mutation-distances to the antigen of more than $`12`$ seem to be unphysiological . Only a cutout of the shape space in the neighborhood of the antigen (which is fixed at $`\mathrm{\Phi }_0`$ in the shape space) is considered: The range of each shape space dimension is $`0\mathrm{\Phi }_i10`$. B cells that leave this shape space area are considered to have vanishing affinity to the antigen. The positions of the seeder cells are listed in Tab. 2. As we use a stochastic model, the GC reactions are repeated with different initializations of the random number generator. Each result is presented by showing $`13`$ GC reactions, the average value, and one standard deviation.

## 3 Results

We find a clear correlation of the total number of output cells produced during the GC reactions and the mutation-distance of the seeder cells that initiate the GC reaction (see Fig. 2). If the seeder cells are already high affinity B cells, the GC reaction is rather stable and produces a large amount of high affinity output cells. With increasing mutation-distance the number of output cells is reduced in a sigmoidal-like functional dependence. For mutation-distances of more than $`6`$ the outcome of the GC reaction becomes critical and for more than $`9`$ the GC reaction is not able to generate a relevant amount of output cells.

The GC reaction breaks down in a relatively narrow range of the seeder cell mutation-distance to the antigen. This points towards an all-or-none development of GC reactions which depends on the seeder cell quality (at least concerning the output quantity). This interpretation is verified by counting the numbers of GC reactions that lead to the production of different numbers of output cells (see Fig. 3). This analysis is based on the same number of GC reactions for each type of seeder cells. Indeed, GCs preferentially produce large or nearly vanishing amounts of output cells. Small or moderate amounts are rare events.

Note, that these results remain basically unaltered considering seeder cells with mutation paths to the antigen of moderate (instead of maximum) probabilistic weight (see Tab. 2, Fig. 2 and 3). However, the break-down of the GC reaction is slightly shifted to seeder cells of higher quality.

The output quality is correlated to the seeder cell quality in a very analogous way (see Fig. 4). The GC reaction produces a stable and large fraction of high affinity output cells for high affinity seeder cells. For seeder cells with larger mutation-distances the fraction of high affinity output decreases. The correlation becomes instable for very large mutation-distances. This is due to the small number of produced output cells in this regime. Large statistical fluctuations obscure the correlation of output and seeder cell quality. However, plotting the number of GC reaction leading to some fraction of high affinity output cells (see Fig. 5) shows that the all-or-none behavior of GCs is clearly found for the output quality as well. Note, that this result is again not affected by switching to seeder cells with mutation paths to the antigen of moderate (instead of maximum) probabilistic weight.

An interesting observation of the GC reactions for various types of seeder cells is that the dark zone is established in all GC reaction in a very similar way – independently of the seeder cell quality. The dark zone is always depleted around day $`8`$ of the reaction (data not shown). This implies that the initial phase of the GC reaction is not influenced by the seeder cell affinity to the antigen. This initial phase includes the phase of monoclonal expansion and even the first days of the phase of somatic hypermutation and selection (what we have called the primary optimization phase ). It is only afterwards when the selection of centrocytes becomes the dominant process of the GC reaction, that the seeder cell affinity to the antigen is a critical parameter. We observe a huge sensitivity of the total GC life time on the seeder cell quality, which varies from $`12`$ to $`27`$ days. This implies that if the B cell quality is not substantially enhanced during the first phase of the reaction, the GC B cell population rapidly dies out within $`4`$ to $`5`$ days – a behavior that has been observed in experiment as well .

It is important to verify that the affinity maturation process was, indeed, unsuccessful in rapidly dying GCs. This is best shown by considering the number of recycled centroblasts (i.e. re-proliferating centroblasts that have been positively selected at least once) after $`10`$ days (a representative time point) of the reaction. If a centrocyte is positively selected it will restart to proliferate with a high probability  and therefore the number of recycled cells is a good measure for the success of the affinity maturation process. For seeder cells with a mutation-distance of more than $`8`$ this value remains below $`20`$ and reaches $`150`$ for seeder cells with mutation-distance of $`3`$.

Let us have a deeper look at the regime where the GC reaction becomes critical. This is the case for seeder cells with mutation-distances from $`6`$ and $`7`$. The same analysis as performed in Figs. 3 and 5 is repeated taking into account the corresponding GC simulations only. The result is shown in Fig. 6. The all-or-none behavior seen before has disappeared. The distribution of GCs with respect to the fraction of high affinity output cells is peaked at $`65\%`$ and flat for smaller fractions. In the case of the number of output cells the distribution of GC reactions is clearly peaked at about 50 cells in a Gaussian-like shape. This is the result expected if the all-or-none behavior of GCs is indeed correlated to the seeder cell quality. If, considering a small window of mutation-distances for the seeder cells, the all-or-none behavior still appeared as it did in Figs. 3 and 5, the reason for the all-or-none behavior would have to be looked for somewhere else in the GC dynamics.

## 4 Discussion

The present analysis has revealed that the seeder cell quality (measured as the mutation-distance to the optimal clone) determines the fate of the GC reaction and is responsible for the all-or-none development of GC reactions as observed in experiments. The efficiency of the GC reaction has been measured quantitatively and qualitatively considering number and affinity of the produced output cells. This result implies that the repertoire of B cells in an organism determines whether the immune system may respond with GC reactions in an efficient way or not. If the repertoire is reduced under some threshold, the probability of a successful GC reaction may become critically small for some specific antigens.

However, there are some model assumptions that should be revisited. At first, the model is two-dimensional. In three dimensions the number of B cells in a fully developed GC is substantially larger. Therefore, we expect a shift of the curves in Fig. 2 and 4 to higher mutation-distances of the seeder cells. This shift is nontrivial and difficult to be estimated quantitatively because the translation from two to three dimension may also change the interaction frequency of centrocytes and FDCs.

In view of the shape space concept a possible role of key mutations may affect the results even when the key mutations are in a cold spot as in the case of NP . However, one may expect that the existence of key mutations would even intensify the all-or-none behavior of GCs because the B cells are confronted with an additional critical process, namely to find the key mutation or not. Therefore, one may suspect that the sigmoidal function in Fig. 2 may even become steeper. From this perspective the present results show that the existence of key mutations is not necessary to understand the all-or-none development of GC reactions . The latter follows already within a quasi-continuous affinity maturation concept.

The analysis is based on the assumption that the seeder cells are homogeneously distributed with respect to the mutation-distances to the antigen. This assumption is not essential for the outcome provided that not all activated B cells that give rise to a GC reaction always have the same mutation-distance to the antigen (see Fig. 6). It is reasonable to expect that the mutation-distance of activated B cells is basically a random parameter. Indeed, it has recently been found for T cell-dependent primary immune responses to NP, that high and low affinity B lymphocytes are equally suitable to induce GC reactions . The affinity of B cells to the antigen differed by a factor of 40 in those measurements. Therefore, we expect that deviations from our assumption of homogeneity will remain sufficiently small not to affect our interpretation of the results.

An important outcome of the analysis is that the seeder cell quality does not affect the initial phase of the GC reaction (compare also ). One may have suspected that the GC reaction simply doesn’t start if the seeder cells have subthreshold affinity to the antigen  (it is assumed that the affinity is sufficient to activate some naive B cells). This possibility is not excluded by the present analysis, as there may be some process, which is not included in the model, but which is crucial for the initiation of the GC reaction. However, our computer experiment shows that the all-or-none behavior can be understood on the level of the efficiency of GC reactions without postulating an additional process during the initiation of the GC reaction.

This result is in accordance with observations of T cell independent GC reactions . Here, it has been found that GCs develop normally in absence of T cell help. Even the dark and light zones are established . As T cells are believed to be of great importance for the selection process (especially to inhibit apoptosis of centrocytes ) the author claimed that the initial phase is basically independent of T cell help. In the present context, this implies that provided that B cells have been activated by the antigen, they may induce a GC reaction independently of their affinity to the antigen (note that assuming a preselection process of un-mutated B cells , this statement applies to positively pre-selected B cells). In a later phase of the GC reaction the T cell help is crucial for centrocytes to survive and, consequently, the already well established GCs are rapidly depleted without T cells . The same happens in the present simulations but for a different reason: While in T cell independent GC reactions the selection process cannot start at all, bad quality seeder cells may still allow for a weak selection process, which is, however, insufficient to achieve a successful affinity maturation process. In both scenarios – in T cell independent GCs and for GCs initiated by low quality seeder cells – the GCs die out because of a failing selection procedure.

The differences of GC responses to two similar (but non-identical) peptides (PS1CT3 and G32CT3) have been examined . The authors found that the intensity of a GC reaction is dependent on the used peptide and, consequently, dependent on the quality of the initially activated B cells. The results found here support these findings, and, in addition, specify the character of this dependence: The seeder cell quality induces an all-or-none behavior of the GC reaction intensity.


text:
# New Experimental Constraints on Non-Newtonian Forces below 100 𝜇m

## Abstract

We have searched for large deviations from Newtonian gravity by means of a finite-frequency microcantilever-based experiment. Our data eliminate from consideration mechanisms of deviation that posit strengths $`10^4`$ times Newtonian gravity at length scales of $`20`$ $`\mu `$m. This measurement is 3 orders of magnitude more sensitive than others that provide constraints at similar length scales.

Physics beyond the standard model will be explored at the Large Hadron Collider in a few years, but gravitational measurements at small length scales are capable of investigating this region now. These measurements become particularly important in light of recent theoretical developments. Moduli, massive scalar particles that could mediate forces with strengths possibly $`10^6`$ times Newtonian gravity, could have Compton wavelengths in the range of $`10`$ $`\mu `$m–$`1`$ mm. These particles would manifest themselves as strong gravitylike forces. The hierarchy problem, the seeming disparity ($`16`$ orders of magnitude) between the standard model energy scale and the Planck scale, suggests other reasons to study this regime. This problem can be recast by asking why gravity is so weak compared to the other known forces in nature. Recent theoretical work add  suggests that extra spatial dimensions, possibly as large as $`1`$ mm and accessible only to gravitons, may supply a solution. Because of spreading into the extra dimensions, gravity would be diluted compared to the other forces, and therefore would seem quite weak. As gravity has not been well tested below the centimeter scale until recently adel ; review , there are many reasons to investigate gravitylike forces at length scales of $`100`$ $`\mu `$m.

Exotic massive scalar particles would modify the standard Newtonian potential with an additional Yukawa-type term scalar . A similar correction would exist for observation of gravity at length scales roughly the size of any compactified extra dimensions add . This leads to the following equation for the gravitational potential in the presence of such non-Newtonian effects:

$$V=G\frac{m_1m_2}{r}(1+\alpha e^{r/\lambda })$$
(1)
Here $`G`$ is Newton’s constant, $`m_1`$ and $`m_2`$ are the masses, $`r`$ is the center of mass separation, and $`\alpha `$ is the strength (relative to Newtonian) of any new effect with a length scale of $`\lambda `$. Non-Newtonian effects are typically parameterized in terms of $`\alpha `$ and $`\lambda `$.

In this Letter we present data obtained for attracting-mass surface separations down to $`25`$ $`\mu `$m, enabling investigation of interaction scales below $`10`$ $`\mu `$m.

The experiment is reminiscent of torsion pendulum gravity measurements, though an ac drive force is employed, and the masses are much smaller. In place of the typical torsional fiber force sensor, a microcantilever was used. A mass attached to this cantilever was subjected to a time-varying force of a gravitational type, and the deflection of the cantilever (measured using fiber-optic interferometry rugar ) provided a measure of this force. The experiment was performed at low temperature and in cryogenic vacuum to exploit the high quality factors attainable in microcantilevers under these conditions.

The force sensors for these experiments were single-crystal silicon diving-board-shaped oscillators $`250`$ $`\mu `$m in length, $`50`$ $`\mu `$m in width, and $`0.335`$ $`\mu `$m thick. The cantilevers were fabricated using standard micromachining techniques atto  and have quality factors in the $`10^310^5`$ range (in vacuum). Spring constants of these sensors were $`5.0`$$`5.5`$ mN/m, and the resonant frequency of the lowest flexural mode of the cantilever was shifted from $`7000`$ Hz to $`300`$ Hz upon mass attachment.

Microcantilever sensitivity is usually limited by thermal noise, i.e., Brownian motion of the beam. This noise limit, analogous to the Johnson noise in a resistor, can be used to obtain the minimum detectable force of a cantilever: $`F_{min}=\sqrt{4k_BTkB/\omega _0Q}`$. Here $`k_B`$ is the Boltzmann constant, $`T`$ is the temperature, $`k`$ is the spring constant, $`\omega _0`$ is the resonant frequency (rad/s), $`Q`$ is the quality factor, and $`B`$ is the measurement bandwidth. The typical force sensitivity for these (thermal noise limited) measurements was $`1\times 10^{16}`$ N/$`\sqrt{\mathrm{Hz}}`$. With averaging times of a few hours, the present devices could reach an ultimate sensitivity below $`1\times 10^{17}`$ N ($`200`$ times the expected Newtonian gravitational force).

The force measured was that between a mass placed on the end of the cantilever (the test mass) and a larger mass oscillated a small distance away (the drive mass). The test masses were gold rectangular prisms, $`50\times 50\times 30`$ $`\mu `$m3 in size and $`1.4`$ $`\mu `$g in mass, and were attached to the cantilevers with a thin (typically $`1`$ $`\mu `$m) layer of epoxy. The drive mass was constructed of ten bars: five gold bars alternating with five silicon bars. Each bar was $`100`$ $`\mu \mathrm{m}\times 100`$ $`\mu \mathrm{m}\times 1`$ mm long. (See Fig. 1 for schematic representation.) This construction, when oscillated in the plane of the drive mass and perpendicular to the direction of the long dimension of the bars, is expected to gravitationally excite the test mass at a harmonic of the oscillation frequency. This shift of the gravitational signal frequency with respect to that of the drive mass oscillation helped prevent spurious excitation of the cantilever due to unwanted vibration at the signal frequency.

To test the system and precisely align the drive mass with the test mass on the cantilever, as well as to provide an in situ equivalent measurement, a magnetic analog of the gravitational experiment was used. The gold bars in the drive mass were connected electrically at alternating ends (see Fig. 1) to form a meander. When a dc electric current is driven through the path defined by the gold bars and their interconnects, a space-varying magnetic field (with a spatial periodicity half that of the gravitational response) is created above the drive mass. This field couples to the test mass through an evaporated layer of nickel on the test-mass surface closest to the drive mass. The oscillation of the drive mass creates a time-varying magnetic field at the location of this layer. A signal proportional to current through the meander verifies the system’s functionality.

To fabricate drive and test masses, bulk silicon was etched to form a mold into which gold was evaporated. After evaporation, the mold was ground and polished. In the case of the test masses, the silicon was then dissolved to release the gold prisms (for fabrication details, see our\_rsi ).

The oscillation of the drive mass beneath the test mass was facilitated by attachment of the drive mass to the free end of a clamped piezoelectric bimorph actuator. The actuator was $`1`$ mm thick and $`4`$ mm wide and had a free length of $`41`$ mm.

The actuator’s motion was characterized ex situ using a laser beam-bounce method. In the experimental probe, the actuator’s motion was measured capacitively using a calibration obtained from this initial characterization. The actuator’s drive frequency was tuned to one-third of the cantilever resonant frequency by means of magnetic excitation. The amplitude of actuator motion was 98$`\pm `$$`\mu `$m. Finite element calculations were performed to determine the expected Newtonian and non-Newtonian response of the test mass.

The nonlinearity in the piezoelectric actuator caused a small amount ($`4\%`$) of actuator vibration at the cantilever resonance (the third harmonic of the actuator frequency) and could lead to a spurious signal. Therefore, two mass-spring vibration isolation stages separated the cantilever mount from the actuator mount. Each stage had a resonant frequency of $`2`$ Hz for both vertical and horizontal motion, and thus together the stages provided an attenuation of $`10^8`$ at $`300`$ Hz between the mounts.

A stiff metallized shield was placed between the drive and test masses to minimize the effects of electrostatic and Casimir forces by preventing direct ac coupling between the masses. The shield was a $`2.8`$ mm by $`5.2`$ mm silicon-nitride membrane $`3`$ $`\mu `$m thick onto both sides of which was evaporated $`100`$ nm of gold. Supported by a 1 cm2 silicon wafer die, the shield was attached to the cantilever wafer $`15`$ $`\mu `$m from the cantilever.

The vertical $`(Z)`$ and in-plane $`(X,Y)`$ relative displacements of the test and drive masses were determined by means of a capacitive position sensor (CPS) similar to that described in field\_cps . Two simpler capacitors, in conjunction with the CPS, provided measurement of relative tilt. After cooling, the position was adjusted until the original capacitive readings were regained, using a three-axis translation stage. Because of the drive mass geometry, the position in $`X`$ (the direction along the drive mass bars) need only be determined to an accuracy of $`100`$ $`\mu `$m. In the $`Y`$ direction, greater precision was required as the gravitational response was expected to vary by $`100`$% over 50 $`\mu `$m in $`Y`$. The CPS enabled determination of relative position in this direction to $`\pm 2`$ $`\mu `$m. The tilt was adjusted by differential heating of the upper vibration isolation stage support springs.

After lateral realignment, an absolute value of the vertical separation between the masses was determined by means of direct mechanical contact (sensed by the cantilever) between the drive mass and the shield. Subsequent motion away from the shield was measured with the CPS, and the separation between the masses was thus known to an accuracy of $`2`$ $`\mu `$m.

The cantilever, actuator, and vibration isolation stages were suspended in a vacuum can at the end of a probe that was inserted into a liquid helium research dewar. An exchange gas space separated the inner vacuum can from the liquid helium in order to reduce external vibrations due to helium boiloff. The entire system was supported by $`2`$ m long $`2`$ Hz springs for additional isolation. The experiments were performed at $`T=9`$$`11`$ K, and at gas pressures of less than $`10^4`$ torr (room temperature reading) with further cryopumping while cold.

The sinusoidal actuator voltage and the cantilever displacement signal were sampled and stored simultaneously at $`10`$ kHz using an analog-to-digital converter. The displacement data stream was then sorted into bins using the phase of the actuator drive voltage and averaged by actuator phase. As the actuator phase is a measure of the lateral position of the drive mass, this technique is similar to a lock-in technique in that signals that are not phase coherent with the actuator drive signal are rejected.

The third harmonic of the Fourier transform of the averaged data is the motion of the cantilever that is phase correlated with the third harmonic of the drive mass motion, i.e., the gravitational-like response of the test mass \[see Fig. 2(a)\]. Using the spring constant and the measured quality factor of the cantilever, the force on the test mass due to the drive mass can be determined.

Data were collected as a function of drive mass meander current and $`Y`$ position in order to align the drive mass to the test mass. Figure 2(b) shows data for one period of the drive mass actuator motion as a function of drive mass meander current. The phase of the magnetic response changes by $`\pi `$ as the current is changed from negative to positive. Also present is a small asymmetry in the magnitude of the magnetic signal for positive and negative values of the nominally same current magnitude. The observed zero-current signal adds phasewise to the signal for a particular current, producing this asymmetry.

The third harmonic of two sets of data acquired over $`1000`$ sec and averaged as a function of measurement time is plotted in Fig. 3 ($`Q=1180`$ for these data). Also plotted is a line representing the theoretical thermal noise of the cantilever. Data taken while the actuator was oscillating far from the test mass match the theoretical curve within error footnote1  and imply that the dominant experimental noise is thermal in nature. (Note that the theoretical curve depends on experimental parameters.) Data taken with the drive mass oscillating and in close proximity to the test mass show a signal that clearly is not diminished with averaging and is patently above the thermal noise. According to finite-element calculations yielding an expected gravitational response of 4.0$`\times 10^{20}`$ N, this anomalous signal is too large (by a factor of $`2000`$) to be of Newtonian origin. Within experimental error, the magnitude of the signal (relative to the actuator drive) as a function of $`Y`$-direction offset is not consistent with a mass-dependent force.

These data lead to the conclusion that the experiment’s force resolution is presently limited by an environmental effect, most likely an electrostatic interaction between the shield and the cantilever as discussed below. While the shield prevents direct coupling of this sort between the masses, a secondary effect due to shield motion could spuriously drive the cantilever. Other possible non-mass-dependent effects are too small to appear at this level. For instance, randomly quenched magnetic impurities in the gold of the drive mass produce a moment that would result in a force below the level of Newtonian gravity, and coordinated pressure variations in the ballistic, residual gas in the probe are similarly small.

Shield motion was characterized (using the fiber interferometer) under standard experimental conditions. The response was not sensitive to the current through the drive mass meander or relative mass position. With the actuator swinging its full amplitude, the shield motion at the first harmonic was $`10`$ pm. The signal at the third harmonic was an order of magnitude less.

Motion of the membrane could drive the cantilever by creating an ac Casimir force. An upper bound can be calculated assuming the cantilever is a conductor and using the Casimir result casimirp  for two conducting planes (force per unit area), $`F_A=\pi ^2\mathrm{}c/240z^4`$. For $`10`$ pm of motion of the shield, the oscillatory force is on the order of $`1\times 10^{21}`$ N, much less than the thermal noise of the cantilever for experimental measurement times.

The motion of the membrane could also induce cantilever motion if there were a significant potential difference between the metallized top of the shield and the cantilever. Making a worst-case assumption of two parallel plates with an area the size of the cantilever, a voltage difference $`0.3`$ V would produce a force on the cantilever of $`1\times 10^{16}`$ N. Potentials of this size have been observed in other experiments stipe , and we cannot yet exclude the possibility that such a force currently limits our experimental sensitivity. Shield motion was most likely due to inhomogeneities in the drive mass surface caused by the use of electrically dissimilar materials. The design of upcoming experiments contains insulating and conducting layers above the drive mass bars; the conducting layer will shield the effect of any drive mass surface variations.

The existence of a spurious force of average magnitude $`8.9\times 10^{17}`$ N allows calculation of an upper limit on the magnitude of any gravitylike force. Figure 4 displays the region of non-Newtonian interaction strength versus interaction distance parameter space that can be excluded by these results footnote2 , as well as other current experimental bounds. For a length scale of $`15`$ $`\mu `$m, the present data improve constraints on the strength of non-Newtonian effects by $`3`$ orders of magnitude. The region eliminated is an important one, as it contains previously unexplored space for scalar moduli particles.

We have presented results from a search for gravitylike deviations from Newtonian theory at distance scales below $`100`$ $`\mu `$m using a microcantilever approach, with masses of size on the order of these scales. Our data offer a new limit on non-Newtonian effects in the range of $`10`$ $`\mu `$m, constraining theoretical predictions of moduli.

We thank Thomas Kenny and Savas Dimopoulos for useful discussions. J.C. and D.M.W. thank DoD for support. S.J.S. thanks NSF for support. This work made use of the National Nanofabrication Users Network facilities supported by the National Science Foundation under Award No. ECS-9731294.

Note added.—After submission of this manuscript, results for larger length scales were published from a similar experiment at the University of Colorado pricenature .