text:
# Kummer structures on a K3 surface - An old question of T. Shioda

## 1. FM partners and Kummer structures

In this section, we shall prove assertion 1 of Main Theorem.

Proof of assertion 1. Let $`A`$ and $`B`$ be abelian surfaces. By \[Ni1\] (see also \[Mo, Lemma 3.1\]), the canonical rational map $`\pi _A:A\text{-}\text{-}\text{-}>\text{Km}A`$ gives the Hodge isometry

$$\pi _{A,}:(T(A)(2),𝐂\omega _A)(T(\text{Km}A),𝐂\omega _{\text{Km}A}),$$

and likewise for $`B`$ and $`\text{Km}B`$. Here, for a lattice $`L:=(L,(,))`$ and for a non-zero integer $`m`$, we define the lattice $`L(m)`$ by $`L(m):=(L,m(,))`$.

Therefore, one has

$$(T(B),𝐂\omega _B)(T(A),𝐂\omega _A)(T(\text{Km}B),𝐂\omega _{\text{Km}B})(T(\text{Km}A),𝐂\omega _{\text{Km}A}).$$

Hence, by the characterization of the FM partners of K3 and abelian surfaces, one obtains

$$BFM(A)\text{Km}BFM(\text{Km}A).$$

On the other hand, we know that $`FM(\text{Km}A)=\{\text{Km}A\}`$ by \[Mu2\] or by the Counting Formula \[HLOY, §3\]. Then we have

$$\text{Km}BFM(\text{Km}A)\text{Km}B\text{Km}A,$$

and therefore

$$BFM(A)B𝒦(X).$$

Since $`|FM(A)|<\mathrm{}`$ by \[BM, Proposition 5.3\] and $`𝒦(X)=FM(A)`$, we have now $`|𝒦(A)|<\mathrm{}`$. See also Section 2 for a more explicit estimate of $`|FM(A)|`$. As it is remarked in Introduction, $`\widehat{A}FM(A)`$. Thus, by $`𝒦(X)=FM(A)`$, we have $`\widehat{A}𝒦(X)`$ as well. ∎

## 2. Kummer structures in larger Picard numbers

In this section, we shall prove assertion 2 of Main Theorem. For this purpose, we shall first give an effective estimate of $`|FM(A)|`$. This is based on the argument of \[HLOY, §3\] and the following Theorem due to Shioda \[Sh2\] (See also \[Mo\]):

###### Theorem 2.1.

1. Let $`A`$ and $`B`$ be complex $`2`$-tori. Then

$$(H^2(B,𝐙),𝐂\omega _B)(H^2(A,𝐙),𝐂\omega _A)BA\mathrm{or}B\widehat{A}.$$
2. For any given weight two Hodge structure $`(U^3,𝐂\omega )`$, there is a marked complex $`2`$-torus $`(A,\tau _A)`$ such that $`\tau _A`$ gives a Hodge isometry

$$\tau _A:(H^2(A,𝐙),𝐂\omega _A)(U^3,𝐂\omega ).$$

For the second statement, we recall that $`U`$ is the even unimodular hyperbolic lattice of rank $`2`$ and the second cohomology lattice $`H^2(A,𝐙)`$ of a complex $`2`$-torus $`A`$ is always isometric to the lattice $`U^3`$.

Let $`A`$ be an abelian surface. Denote by $`G=O_{Hodge}(T(A),𝐂\omega _A)`$ the group of Hodge isometries of $`(T(A),𝐂\omega _A)`$. We call two primitive embeddings $`\iota :T(A)U^3`$, $`\iota ^{}:T(A)U^3`$ $`G`$-equivalent if there exist $`\mathrm{\Phi }O(U^3)`$ and $`gG`$ such that $`\iota ^{}g=\mathrm{\Phi }\iota `$, i.e. the following diagram commutes;

$$\begin{array}{ccc}\iota :T(A)& & U^3\\ g& & \mathrm{\Phi }& & \\ \iota ^{}:T(A)& & U^3\end{array}$$

We denote by $`𝒫^{G\text{-eq}}(T(A),U^3)`$ the set of $`G`$-equivalence classes of primitive embeddings of $`T(A)`$ into $`U^3`$:

$$𝒫^{G\text{-eq}}(T(A),U^3):=\{\text{primitive embedding }\iota :T(A)U^3\}/G\text{-equiv.}.$$

Claim. There is a surjection,

$$\xi :FM(A)𝒫^{G\text{-eq}}(T(A),U^3),B[\iota _B],$$

such that $`\xi ^1([\iota _B])=\{B,\widehat{B}\}`$.

###### Proof.

Let $`BFM(A)`$. Choose a marking $`\tau _B:H^2(B,𝐙)U^3`$ and a Hodge isometry $`\gamma :(T(A),𝐂\omega _A)(T(B),𝐂\omega _B)`$. Define $`\xi :[B][\iota _B]`$, where $`\iota _B=\tau _B\gamma `$. This is a well-defined surjection with $`\xi ^1(\xi (A))=\{A,\widehat{A}\}`$. The proof of this fact is exactly the same as the argument for Theorem 3.3 of \[HLOY\], except that Theorem 2.1 here plays the corresponding role of the Torelli Theorem for K3 surfaces there. Note that $`\xi `$ may fail to be injective because $`B`$ and $`\widehat{B}`$ need not be isomorphic in general. ∎

Recall Theorem 1.4 of \[HLOY\]:

(2.1)
$$|𝒫^{G\text{-eq}}(T(A),U^3)|=\underset{j=1}{\overset{m}{}}|O(S_j)O(A_{S_j})/G|.$$

Here $`𝒢(NS(A)):=\{S_1,S_2,\mathrm{},S_m\}(S_1NS(A))`$ is the genus of the Néron-Severi lattice $`NS(A)`$ and $`O(A_{S_i})`$ is the orthogonal group of the discriminant group $`A_{S_i}=S_i^{}/S_i`$ with respect to the natural $`𝐐/2𝐙`$-valued quadratic form $`q_{S_i}`$ and the set $`O(S_j)O(A_{S_j})/G`$ is the orbit space of the natural action of $`O(S_j)\times G`$ on $`A_{S_j}`$. Note that $`NS(A)`$ is even and hyperbolic. We also recall that the genus of $`NS(A)`$ is nothing but the set of isomorphism classes of even lattices $`S`$ such that $`(NS(A),q_{NS(A)})(S,q_S)`$, by \[Ni2, Corollary 1.9.4\].

We can now prove assertion 2 of Main Theorem.

If $`\text{det}NS(A)`$ is square free, it follows that $`A_{NS(A)}`$ is cyclic, i.e. $`l(A_{NS(A)})=1`$. Here $`l(M)`$ denote the minimal number of generators of a finite abelian group $`M`$. Hence, if, in addition, $`\rho (A)=3`$, then one has

$$\text{rank}NS(A)=32+l(A_{NS(A)}).$$

Similarly, if $`\rho (A)=4`$, then one has

$$\text{rank}NS(A)=42+l(A_{NS(A)}).$$

Indeed, if $`\rho (A)=4`$, then $`\text{rank}T(A)=2`$. Since $`(A_{NS(A)},q_{NS(A)})(A_{T(A)},q_{T(A)})`$, we have $`l(A_{NS(A)})=l(A_{T(A)})\text{rank}T(A)=2`$.

Hence, by \[Ni2, Theorem 1.14.2\], one has that $`𝒢(NS(A))=\{NS(A)\}`$ and that the natural map $`O(NS(A))O(A_{NS(A)})`$ is surjective in each case. Now it follows from the equation (2.1) that $`|𝒫^{G\text{-eq}}(T(A),U^3)|=1`$. Hence, by the preceding claim, we have $`FM(A)=\{A,\widehat{A}\}`$. Combining this with assertion 1 of Main Theorem, we obtain the desired equality $`𝒦(X)=\{A,\widehat{A}\}`$. ∎

Remark. In \[SM\], it is shown that an abelian surface $`A`$ with $`\rho (A)=4`$ is necessarily a direct product of two elliptic curves, namely, $`AE\times F`$. Since $`\widehat{E}E`$ and $`\widehat{F}F`$, it follows that $`\widehat{A}A`$ if $`AE\times F`$. Combining this with assertion 2, one has $`𝒦(X)=\{A\}`$ if $`\rho (A)=4`$. This gives an alternative explanation for the result \[SM, Theorem 5.1\] quoted in Introduction. ∎

## 3. K3 surfaces with many Kummer structures

In this section, we shall show assertion 3 of Main Theorem.

We begin by the following result first found in \[Og2\]:

###### Theorem 3.1.

For any given natural number $`N`$, there is an even hyperbolic lattice $`S`$ such that $`\text{rank}S=2`$, $`\text{det}S=pq`$ and $`|𝒢(S)|N`$. Here $`𝒢(S)`$ is the genus of the lattice $`S`$, and $`p`$ and $`q`$ are some distinct prime numbers (including $`1`$).

Remark. Such phenomena seem quite rare and are perhaps impossible if $`\text{rank}S3`$ (cf. \[Cs, Chapters 9 and 10\] and \[CL, §8 - §10\]). In this sense, as it is pointed out to us by Professor B. Gross, this Theorem is a sort of miracle which happens probably only in rank $`2`$. ∎

###### Proof.

Here we shall sketch the argument. Our argument requires some basic properties about class groups. These are found in the book \[Za\]. (See also \[HLOY, Section 4\] for summary.) Let $`P`$ be the set of prime numbers including $`1`$. By \[Iw, Main Theorem\], the set

$$𝒩:=\{(n,p,q)𝐍\times P^2|4n^2+1=pq\}$$

is an infinite set. Note that $`pq1\text{mod}\mathrm{\hspace{0.17em}4}`$ and $`pq`$ for $`(n,p,q)𝒩`$. Therefore $`pq`$ is a so-called fundamental discriminant. Let $`H(pq)`$ be the narrow class group of (the ring of integers $`O(pq)`$ of) the real quadratic field $`𝐐(\sqrt{pq})`$. We denote by $`h(pq)`$ the class number of $`𝐐(\sqrt{pq})`$, i.e. $`|H(pq)|`$. Recall that there is a natural bijective correspondence between $`H(pq)`$ and the set $`(pq)`$ of the proper, i.e. $`\text{SL}(2,𝐙)`$, isomorhism classes of even hyperbolic rank $`2`$ lattices $`S`$ with $`detS=pq`$. Recall also that $`(pq)`$ is decomposed into at most two sets of the same cardinality:

$$(pq)=\stackrel{~}{𝒢}(S_1)\stackrel{~}{𝒢}(S_2).$$

Here $`\stackrel{~}{𝒢}(S_i)`$ is the set of proper, i.e. $`\text{SL}(2,𝐙)`$, isomorphism classes of lattices in the same genus as $`S_i`$. Combining these two facts with $`\text{GL}(2,𝐙)/\text{SL}(2,𝐙)𝐙/2`$, we have

$$|𝒢(S_1)|\frac{|\stackrel{~}{𝒢}(S_1)|}{2}\frac{|(pq)|}{4}=\frac{h(pq)}{4}.$$

Therefore, it is enough to show that $`h(pq)\mathrm{}`$ when $`pq\mathrm{}`$. We note that there is a sequence $`\{(n_k,p_k,q_k)\}_{k=1}^{\mathrm{}}𝒩`$ such that $`p_kq_k\mathrm{}`$ when $`k\mathrm{}`$. This is because $`|𝒩|=\mathrm{}`$.

Let us show that $`h(pq)\mathrm{}`$ when $`pq\mathrm{}`$. Since $`(2n+\sqrt{pq})(2n+\sqrt{pq})=1`$ by $`4n^2+1=pq`$, the element $`2n+\sqrt{pq}>1`$ is a unit of $`O(pq)`$. Therefore, the fundamental unit $`ϵ(pq)`$ of $`O(pq)`$ satisfies $`1<ϵ(pq)<pq`$. Now, combining this with the Siegel-Brauer formula

$$\underset{pq\mathrm{}}{lim}\frac{\mathrm{log}(h(pq)\mathrm{log}ϵ(pq))}{\mathrm{log}pq}=\frac{1}{2},$$

one has $`h(pq)\mathrm{}`$ when $`pq\mathrm{}`$. ∎

Let us return to the proof of assertion 3 of Main Theorem.

Let $`S`$ be the lattice in Theorem (4.1) and take $`N`$ non-isomorphic elements

$$S_1:=S,S_2,\mathrm{},S_N𝒢(S).$$

Then $`S_i`$ are even hyperbolic lattices of $`\text{rank}S_i=2`$ and satisfy $`(A_{S_i},q_{S_i})(A_S,q_S)`$.

Let us choose a primitive embedding

$$\phi _i:S_iU^3$$

for each $`i`$. For instance, if $`S_i=𝐙v_{i1},v_{i2}`$ with

$$((v_{i1},v_{i2}))=\left(\begin{array}{cc}2a_i& b_i\\ b_i& 2c_i\end{array}\right),$$

then we can define $`\phi _i`$ by $`\phi _i(v_{i1})=e_1+a_if_1`$ and $`\phi _i(v_{i2})=b_if_1+e_2+c_if_2`$, where $`e_j`$ and $`f_j`$ are the standard basis of the $`j`$-th factor $`U`$ of $`U^3`$.

Set $`T_i:=\phi _i(S_i)^{}`$ in $`U^3`$ and $`T:=T_1`$. Then

$$(A_{T_i},q_{T_i})(A_{S_i},q_{S_i})(A_S,q_S)$$

and the signature of $`T_i`$ is $`(2,2)`$ for each $`i`$. Therefore $`T_i𝒢(T)`$.

Since $`pq`$ is square free, we have $`A_TA_S𝐙/pq`$. Thus,

$$\text{rank}T=42+1=2+l(A_T).$$

In addition, $`T`$ is indefinite. Hence $`𝒢(T)=\{T\}`$ by \[Ni2, Theorem 1.14.2\]. Therefore $`T_iT`$ for each $`i`$. Let us choose an isometory

$$\sigma _i:TT_i$$

for each $`i`$. We set $`\sigma _1:=id`$.

Choose a maximal weight two Hodge structure $`(T,𝐂\omega )`$ on $`T`$. Here the term maximal means that $`T`$ satisfies the following property: if $`T^{}T`$ is a primitive sublattice of $`T`$ such that $`𝐂\omega T^{}𝐂`$ then $`T^{}=T`$. This is equivalent to that $`\text{Re}\omega ,\text{Im}\omega ^{}T=\{0\}`$ in $`T𝐑`$. As it is easily shown, there certainly exists a maximal weight two Hodge structure on $`T`$. Note that the Hodge structure $`(T(V),𝐂\omega _V)`$ on $`T(V)`$ of an abelian surface $`V`$ is maximal in this sense.

Set $`\omega _i:=\sigma _i(\omega )`$. Then $`(T_i,𝐂\omega _i)`$ is a maximal weight two Hodge structure on $`T_i`$ such that $`\sigma _i:(T,𝐂\omega )(T_i,𝐂\omega _i)`$. Then, by Theorem 2.1, there is a marked complex $`2`$-torus $`(A_i,\tau _i)`$ such that

$$\tau _i:(H^2(A_i,𝐙),𝐂\omega _{A_i})(U^3,𝐂\omega _i).$$

Then, by the maximality of $`T(A_i)`$ and $`T_i`$, one has the Hodge isometry

$$\tau _i|_{T(A_i)}:(T(A_i),𝐂\omega _{A_i})(T_i,𝐂\omega _i)$$

and whence an isometry of the orthogonal lattices

$$\tau _i|_{NS(A_i)}:NS(A_i)S_i$$

for each $`i`$. Since $`S_i`$ is hyperbolic, the complex $`2`$-tori $`A_i`$ are actually projective, i.e. abelian surfaces. Since $`S_i\simeq ̸S_j(ij)`$ by the choice of $`S_i`$, one has $`NS(A_i)\simeq ̸NS(A_j)`$ ($`ij`$) as well.

Set $`A:=A_1`$. Then, composing the Hodge isometries $`(\tau |_{T(A)}^1)`$, $`\sigma _i^1`$ and $`\tau _i|_{T(A_i)}`$, one has a Hodge isometry

$$(T(A_i),𝐂\omega _{A_i})(T(A),𝐂\omega _A).$$

Therefore $`A_iFM(A)`$ by the characterization of $`FM(A)`$. Thus, $`A_i𝒦(\text{Km}A)`$ by assertion 1 of Main Theorem. Hence, these $`N`$ abelian surfaces $`A_i`$ ($`1iN`$) satisfy the requirement in assertion 3 of Main Theorem. ∎

## Appendix A

In this appendix, we shall show the following:

###### Proposition A.1.

Consider the complex $`2`$-torus $`A=𝐂^2/L`$ defined by a rank $`4`$ discrete lattice $`L`$ in $`𝐂^2`$ which is generated by

$$e_1:=\left(\begin{array}{c}1\\ 0\end{array}\right),e_2:=\left(\begin{array}{c}0\\ 1\end{array}\right),e_3:=\left(\begin{array}{c}\alpha \\ \beta \end{array}\right),e_4:=\left(\begin{array}{c}\gamma \\ \delta \end{array}\right).$$

Assume that $`L`$ is generic in the sense that $`\alpha `$, $`\beta `$, $`\gamma `$ and $`\delta `$ are algebraically independent over $`𝐐`$. Then:

1. $`NS(A)=\{0\}`$ and $`A`$ is not projective.
2. $`\widehat{A}\simeq ̸A`$ and $`\widehat{A}FM(A)`$.
3. $`𝒦(\mathrm{Km}A)=\{A\}`$.

In particular, $`\widehat{A}𝒦(\mathrm{Km}A)`$, i.e. $`\mathrm{Km}\widehat{A}\simeq ̸\mathrm{Km}A`$, while $`\widehat{A}FM(A)`$.

This proposition shows certain differences between the Kummer structures on a non-projective K3 surface and a projective K3 surface.

Proof of assertion 1. Let $`z_1`$ and $`z_2`$ be the coordinates of $`𝐂^2`$ corresponding to the first and the second projections $`𝐂^2𝐂`$. We identify $`L`$ with $`H_1(A,𝐙)`$. Let $`v_i_{i=1}^4`$ be the dual basis of $`H^1(A,𝐙)`$ of the basis $`e_i_{i=1}^4`$ of $`H_1(A,𝐙)`$. Then, by the description of $`e_i`$, we have

$$dz_1=v_1+\alpha v_3+\gamma v_4,dz_2=v_2+\beta v_3+\delta v_4$$

in $`H^1(A,𝐂)`$. From this, we obtain

$$\omega :=dz_1dz_2=v_1v_2+\beta v_1v_3+\delta v_1v_4\alpha v_2v_3\gamma v_2v_4+(\alpha \delta \beta \gamma )v_3v_4.$$

Note that $`NS(A)=\{\eta H^2(A,𝐙)|(\eta ,\omega )=0\}`$. This due to the Lefschetz $`(1,1)`$ Theorem. Let $`\eta NS(A)`$ and write $`\eta `$ as

$$\eta =av_1v_2+bv_1v_3+cv_1v_4+dv_2v_3+ev_2v_4+fv_3v_4.$$

Here the coefficients $`a`$, $`b`$, $`\mathrm{}`$, $`f`$ are integers. By $`(\eta ,\omega )=0`$, one has

$$a(\alpha \delta \beta \gamma )+b\gamma c\alpha +d\delta e\beta +f=0.$$

Since $`\alpha `$, $`\beta `$, $`\gamma `$, $`\delta `$ are algebraically independent over $`𝐐`$, we have then

$$a=b=c=d=e=f=0,\mathrm{i}.\mathrm{e}.,\eta =0.$$

Hence $`NS(A)=\{0\}`$. This also implies that $`A`$ is not projective. ∎

Proof of assertion 2. The second statement $`\widehat{A}FM(A)`$ follows from \[Mu1\], that the Poincaré bundle induces the equivalence $`D(\widehat{A})D(A)`$.

Let us show the first statement. Let $`\widehat{L}`$ be the rank 4 discrete lattice in $`𝐂^2`$ generated by

$$f_1:=\left(\begin{array}{c}1\\ 0\end{array}\right),f_2:=\left(\begin{array}{c}0\\ 1\end{array}\right),f_3:=\left(\begin{array}{c}\alpha \\ \gamma \end{array}\right),f_4:=\left(\begin{array}{c}\beta \\ \delta \end{array}\right).$$

Then, by the formula \[Sh2, (2.13)\], we have $`\widehat{A}=𝐂^2/\widehat{L}`$. Assume that $`A\widehat{A}`$. Then, there is $`B\text{GL}(2,𝐂)`$ such that $`Be_1,Be_2,Be_3,Be_4\widehat{L}`$.

By $`Be_1,Be_2\widehat{L}`$, there are integers $`m_i`$, $`n_i`$ such that

$$Be_1=\left(\begin{array}{c}m_1+m_3\alpha +m_4\beta \\ m_2+m_3\gamma +m_4\delta \end{array}\right),$$

$$Be_2=\left(\begin{array}{c}n_1+n_3\alpha +n_4\beta \\ n_2+n_3\gamma +n_4\delta \end{array}\right).$$

Then, we have

$$B=\left(\begin{array}{cc}m_1+m_3\alpha +m_4\beta & n_1+n_3\alpha +n_4\beta \\ m_2+m_3\gamma +m_4\delta & n_2+n_3\gamma +n_4\delta \end{array}\right).$$

Therefore

$$Be_3=\left(\begin{array}{c}m_1\alpha +m_3\alpha ^2+m_4\alpha \beta +n_1\beta +n_3\alpha \beta +n_4\beta ^2\\ m_2\alpha +m_3\alpha \gamma +m_4\alpha \delta +n_2\beta +n_3\beta \gamma +n_4\beta \delta \end{array}\right).$$

Since $`Be_3\widehat{L}`$, there are integers $`l_i`$ such that

$$\left(\begin{array}{c}m_1\alpha +m_3\alpha ^2+m_4\alpha \beta +n_1\beta +n_3\alpha \beta +n_4\beta ^2\\ m_2\alpha +m_3\alpha \gamma +m_4\alpha \delta +n_2\beta +n_3\beta \gamma +n_4\beta \delta \end{array}\right)=\left(\begin{array}{c}l_1+l_3\alpha +l_4\beta \\ l_2+l_3\gamma +l_4\delta \end{array}\right).$$

Since $`\alpha `$, $`\beta `$, $`\gamma `$, and $`\delta `$ are algebraically independent over $`𝐐`$, this implies that

$$m_2=m_3=m_4=n_2=n_3=n_4=0.$$

However, then $`B`$ is of the form,

$$B=\left(\begin{array}{cc}a& b\\ 0& 0\end{array}\right),$$

and therefore $`\text{det}B=0`$, a contradiction. Hence $`A\simeq ̸\widehat{A}`$. ∎

Proof of assertion 3. Let $`B`$ be a complex $`2`$-torus such that $`\text{Km}B\text{Km}A`$. Let us choose an isomorphism $`\phi :\text{Km}A\text{Km}B`$. Since $`NS(A)=\{0\}`$, there is no effective curves on $`A`$. Therefore $`\text{Km}A`$ contains no effective curves other than the $`16`$ exceptional curves $`C_i`$ arising from the resolution $`\text{Km}AA/1`$. Therefore $`\text{Km}B`$ also contains only 16 $`𝐏^1`$s, say $`D_i`$, by $`\text{Km}B\text{Km}A`$. These $`D_i`$ must then coincide with the $`16`$ exceptional curves arising from the resolution $`\text{Km}BB/1`$ and satisfy $`\phi (_{i=1}^{16}C_i)=_{i=1}^{16}D_i`$. Now recall that $`A`$ can be reconstructed from $`\text{Km}A`$ as follows: first take the (unique) double cover of $`\text{Km}A`$ branched along the divisor $`_{i=1}^{16}C_i`$ and then contract the proper transform of $`_{i=1}^{16}C_i`$. Likewise $`B`$ from $`\text{Km}B`$ and $`_{i=1}^{16}D_i`$ . Since $`\phi (_{i=1}^{16}C_i)=_{i=1}^{16}D_i`$, it follows that $`\phi `$ then lifts to an isomorphism $`\stackrel{~}{\phi }:AB`$. This completes the proof. ∎

Remark. In \[Ni1, main result\], there is a bijective correspondence between the set of Kummer structures $`𝒦(X)`$ on $`X`$ and $`𝒞(X)/\text{Aut}(X)`$, where $`𝒞(X)`$ is the set of the configurations of $`16`$ disjoint $`𝐏^1`$s on $`X`$. However $`𝒞(X)`$ is an infinite set and the natural action of $`\text{Aut}(X)`$ on $`𝒞(X)`$ is not transitive in general, when $`X`$ is a projective K3 surface. It seems much harder to prove our main theorem by such a geometric method using this bijective correspondence.


format: text

text:
# Classical Behaviour After a Phase Transition: I. Classical Order Parameters

## I Introduction

The emergence of classical behaviour from a quantum system is a problem of interest in many branches of physics , from the foundations of quantum mechanics, condensed matter and quantum optics, to quantum computing and quantum field theory. In the last few years, experimental evidence supports our theoretical understanding of the quantum to classical transition and also opens new avenues in the development of essential tools to understand the more hidden mysteries of the quantum world .

The quantum to classical transition, for a point particle, say, involves two different but very related conditions. The first one is that there should be correlations, i.e., the coordinates and momenta of a particle should be correlated in phase space, according to the classical equations of motion. For example, the Wigner function should have a peak at the classical trajectories. The second and equally important condition is the elimination of quantum interference between these classical trajectories, i.e., decoherence. Once decoherence eliminates interference terms (absent in the ‘classical world’), the Wigner function becomes a good candidate for a classical probability distribution. The generalisation from particles to fields is straightforward, in principle, with the same attributes of correlations and decoherence.

The onset of classical behaviour is a natural consequence of a quantum open system, triggered by the interaction between the system of interest and its environment. The coupling strength between system and bath sets the time-scale after which we can consider our system as classical, according to both of the conditions we mentioned before. This temporal scale is usually called the decoherence time $`t_\mathrm{D}`$. After this decoherence time we do not have macroscopic states in coherent superpositions anymore, and a probability distribution can be extracted that evolves by means of a generalised Fokker-Planck equation.

Our concern in this paper will be the quantum to classical transition of a single scalar order-parameter field during continuous transitions, with the simplest double-well potential. This is an idealisation of the phase transitions that are expected to occur at the GUT and EW scales in the standard big-bang cosmological model .

An analysis of phase transitions in quantum field theory that takes the non-equilibrium nature of the dynamics into account from first principles has only recently begun to be addressed. In particular, the naive picture of a classical order parameter (inflaton or Higgs) field $`\varphi `$ rolling down an adiabatic effective potential, that was once a mainstay of cosmological field theory modelling, has been shown to be suspect . Alternatively, the suggestion by Kibble that, while a non-adiabatic approach is essential, causality alone can set saturated bounds on time and distance scales during a transition, has been shown to be only partly true. The onset of classical behaviour is absolutely crucial in the work of Ref. in that the experimental signal for the causal bound is the production of classical defects at a prescribed density. Topological defects are inevitable in most transitions, and they may have played a fundamental role in the formation of large scale structure (strings) . Moreover, superabundance of some topological defects may contradict the observational evidence (magnetic monopoles).

In all the above-mentioned examples there is an order parameter which evolves from the false to the true vacuum of the theory: the Higgs fields in GUT and EW phase transitions, the inflaton field(s) in inflationary models, etc. Although these are quantum scalar fields with vanishing mean value (due to the symmetry of the initial quantum state), the order parameter is usually treated as a classical object. Our aim in this paper is to justify this assumption.

Our approach follows the analysis started by two of us in Ref., where we studied the emergence of classical inhomogeneities from quantum fluctuations for a self-interacting quantum scalar field. We have investigated there the decoherence induced on the long-wavelength field modes by coarse-graining the field modes with wavelength shorter than a critical value, in order to show how the system becomes classical due to the interaction with its environment (in that case composed of the short-wavelength field modes of the same field). For phase transitions the classicality of the order parameter can be analyzed along the same lines by extending the model to accommodate spontaneous symmetry breaking.

This is a difficult problem because, as has been pointed out in the literature, and as we will stress in what follows, non-perturbative and non-Gaussian effects are relevant in the analysis of phase transitions. As a trial run for this analysis we begin by presenting a toy model in which we will study the spread of a wave packet initially centered around the local maximum of a double well potential, paying particular attention to the influence of the environment on the Wigner function and on the reduced density matrix. Only later will we will extend our results to quantum field theory phase transitions.

The paper is organized as follows. In the next section we study the evolution of a wave packet initially centered on the top of a double well potential. We describe the exact numerical evaluation of the evolution of the wave packet. We show that, as the coupling between the system and the environment decreases, the decoherence time increases. Due to the nonlinearities of the potential, when the coupling vanishes there is no classical limit, not even classical correlations. In section III we analyze the full problem in quantum field theory. Here, the situation is somewhat different. Although, again, the decoherence time increases as the coupling to the environment decreases, so does the spinodal time, the time for the order parameter to achieve its classical ground-state values. Thus, however weak the coupling, decoherence can still occur before the field has fallen to its classical level.

It is important to say that the issue of how the system evolves into the classical theory has been addressed in Refs.. For some models, it has been shown that classicality emerges as a consequence of profuse particle creation. A non-perturbative large occupation number of long-wavelength modes produces, on average, a diagonal density matrix. This dephasing effect occurs at late times. Here we will consider a different model in which classicality is an early time event.

Finally, in section IV, we present our final remarks and conclusions.

## II Toy model: double well potential and environment

The simplest field theory that permits a phase transition is that of a single real scalar field $`\varphi `$, with action

$$S[\varphi ]=d^4x\left\{\frac{1}{2}_\mu \varphi ^\mu \varphi +\frac{1}{2}\mu ^2\varphi ^2\frac{\lambda }{4!}\varphi ^4\right\}$$
(1)
with ($`\mu ^2>0`$) $`𝒵_2`$ symmetry breaking. At high temperature the symmetry is restored. As we shall see later, a phase transition induced by a sudden temperature quench can be described by an effective field theory in which there is a change of sign in the mass term of the scalar field

$$S^{\mathrm{eff}}[\varphi ]=d^4x\left[\frac{1}{2}_\mu \varphi ^\mu \varphi \frac{1}{2}m^2(t)\varphi ^2\frac{1}{4!}\lambda \varphi ^4\right],$$
(2)
where $`m^2(t)=M^2>0`$ for $`t<0`$ and $`m^2(t)=\mu ^2`$ for $`t`$ sufficiently positive. This change of sign in $`m^2(t)`$ breaks the global $`𝒵_2`$ symmetry for positive $`t`$.

Understanding this transition, for even such a simple system, is difficult. As a preliminary exercise, we start by considering a toy model composed of a particle, a quantum anharmonic oscillator (the ‘system’), linearly coupled to an environment composed of an infinite set of harmonic oscillators. The total classical action is given by

$`S[x,q_n]`$ $`=`$ $`S_{\mathrm{syst}}[x]+S_{\mathrm{env}}[q_n]+S_{\mathrm{int}}[x,q_n]`$ (3)
$`=`$ $`{\displaystyle _0^t}ds[{\displaystyle \frac{1}{2}}M(\dot{x}^2\mathrm{\Omega }_0^2(t)x^2{\displaystyle \frac{\lambda }{4}}x^4)`$ (4)
$`+`$ $`{\displaystyle \underset{n}{}}{\displaystyle \frac{1}{2}}m_n(\dot{q}_n^2\omega _n^2q_n^2)]{\displaystyle }_nC_nxq_n,`$ (5)
where $`x`$ and $`q_n`$ are the coordinates of the particle and the oscillators, respectively. The quantum anharmonic oscillator is coupled linearly to each oscillator in the bath with strength $`C_n`$. In analogy with $`m^2(t)`$ of (1) we consider the simplest possible case of an instantaneous quench, in which $`\mathrm{\Omega }_0^2(t)=\mathrm{\Omega }_0^2`$ when $`t<0`$ and $`\mathrm{\Omega }_0^2(t)=\mathrm{\Omega }_0^2`$ when $`t>0`$. The unstable particle, coordinate $`x`$, has an initial thermal distribution at temperature $`T`$, for which $`\widehat{x}=0`$. For $`t>0`$ it finds itself on the unstable hump of the potential and falls towards its stable minima at $`\widehat{x}^2=2\mathrm{\Omega }_0^2/\lambda `$.

### A The environment

In the absence of an environment, and considering $`\lambda =0`$, this situation has been discussed in detail by Guth and Pi. In the presence of an environment the relevant objects for analyzing the quantum to classical transition in this model are the reduced density matrix, and the associated Wigner function

$$\rho _\mathrm{r}(x,x^{},t)=𝑑q_n\rho (x,q_n,x^{},q_n,t),$$
(6)
and

$$W_\mathrm{r}(x,p,t)=\frac{1}{2\pi }_{\mathrm{}}^+\mathrm{}𝑑ye^{ipy}\rho _\mathrm{r}(x+\frac{y}{2},x\frac{y}{2},t).$$
(7)
The reduced density matrix satisfies a closed master equation. This has been evaluated by Hu-Paz-Zhang  for the quantum Brownian motion problem with $`\mathrm{\Omega }_0^2`$ positive. Following the same procedure, we can write a corresponding master equation for the unstable particle  by replacing $`\mathrm{\Omega }_0`$ by $`i\mathrm{\Omega }_0`$ in the Hu-Paz-Zhang result

$`i_t\rho _\mathrm{r}(x,x^{},t)`$ $`=`$ $`x|[H,\rho _\mathrm{r}]|x^{}`$ (8)
$``$ $`i\gamma (t)(xx^{})(_x_x^{})\rho _\mathrm{r}(x,x^{},t)`$ (9)
$`+`$ $`f(t)(xx^{})(_x+_x^{})\rho _\mathrm{r}(x,x^{},t)`$ (10)
$``$ $`iD(t)(xx^{})^2\rho _\mathrm{r}(x,x^{},t).`$ (11)
In (11), $`H=H_{\mathrm{syst}}\frac{1}{2}M\stackrel{~}{\mathrm{\Omega }}^2(t)`$ where $`\stackrel{~}{\mathrm{\Omega }}^2(t)`$ renormalizes the natural frequency of the particle, $`\gamma (t)`$ is the dissipation coefficient, and $`D(t)`$ and $`f(t)`$ are the diffusion coefficients, which produce the decoherence effects. They depend on the properties of the environment as

$`\stackrel{~}{\mathrm{\Omega }}^2(t)`$ $`=`$ $`{\displaystyle \frac{2}{M}}{\displaystyle _0^t}𝑑t^{}\mathrm{cosh}(\mathrm{\Omega }_0t^{})\eta (t^{})`$ (12)
$`\gamma (t)`$ $`=`$ $`{\displaystyle \frac{1}{2M\mathrm{\Omega }_0}}{\displaystyle _0^t}𝑑t^{}\mathrm{sinh}(\mathrm{\Omega }_0t^{})\eta (t^{})`$ (13)
$`D(t)`$ $`=`$ $`{\displaystyle _0^t}𝑑t^{}\mathrm{cosh}(\mathrm{\Omega }_0t^{})\nu (t^{})`$ (14)
$`f(t)`$ $`=`$ $`{\displaystyle \frac{1}{M\mathrm{\Omega }_0}}{\displaystyle _0^t}𝑑t^{}\mathrm{sinh}(\mathrm{\Omega }_0t^{})\eta (t^{}),`$ (15)
where $`\eta (t)`$ and $`\nu (t)`$ are the dissipation and noise kernels, respectively,

$`\eta (t)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑\omega I(\omega )\mathrm{sin}\omega t`$ (16)
$`\nu (t)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑\omega I(\omega )\mathrm{coth}{\displaystyle \frac{\beta \omega }{2}}\mathrm{cos}\omega t,`$ (17)
and $`I(\omega )=𝒪(C_n^2)`$ is the spectral density of the environment.

The first term on the RHS of Eq.(11) gives the usual Liouville-like evolution; the term proportional to $`\gamma `$ produces dissipation ($`\gamma `$ is the relaxation coefficient). The term proportional to the diffusion coefficient $`D(t)`$, which is proportional to $`(xx^{})^2`$ and positive definite, gives the main contribution to the decoherence since it produces a diagonalization of the reduced density matrix. Let us write the reduced density matrix as

$$\rho _\mathrm{r}[x,x^{};t]=G[x,x^{},t]\mathrm{exp}\left[(xx^{})^2_0^tD(s)𝑑s\right].$$
(18)
Inserting this expression into the master equation it is easy to to see that the differential equation for $`G[x,x^{},t]`$ contains the usual Liouville-term plus additional contributions proportional to $`D,\gamma `$, and $`f`$. However, none of these additional terms is imaginary with the right sign for diffusion. An approximate solution of Eq.(11) is therefore

$$\rho _\mathrm{r}[x,x^{};t]\rho _\mathrm{r}^\mathrm{u}[x,x^{},t]\mathrm{exp}\left[(xx^{})^2_0^tD(s)𝑑s\right],$$
(19)
where $`\rho _\mathrm{r}^\mathrm{u}`$ takes into account the unitary evolution.

Alternatively, one can derive the following evolution equation for the reduced Wigner function of the system :

$`\dot{W_\mathrm{r}}(x,p,t)`$ $`=`$ $`\{H_{\mathrm{syst}},W_\mathrm{r}\}_{\mathrm{PB}}{\displaystyle \frac{\lambda }{4}}x_{ppp}^3W_\mathrm{r}+2\gamma (t)_p(pW_\mathrm{r})`$ (20)
$`+`$ $`D(t)_{pp}^2W_\mathrm{r}f(t)_{px}^2W_\mathrm{r}.`$ (21)
Let us concentrate on the evolution equation (21). The first term on the right-hand side of Eq.(21) is the Poisson bracket, corresponding to the usual classical evolution. The second term includes the quantum correction (we have set $`\mathrm{}=1`$). The last three terms describe dissipation and diffusion effects due to coupling to the environment. In order to simplify the problem, we consider a high-temperature ohmic ($`I(\omega )\omega `$) environment. In this approximation the coefficients in Eq.(21) become constants: $`\gamma (t)=\gamma _0`$, $`f1/T`$, and $`D=2\gamma _0k_\mathrm{B}T`$. The normal diffusion coefficient $`D`$ is the term responsible for decoherence effects and at high temperatures is much larger than $`\gamma _0`$ and $`f`$. Therefore, in Eq.(21), we may neglect the dissipation and the anomalous diffusion terms in comparison to the normal diffusion. As we saw in (19), it is the diffusion exponential in this equation that enforces the approximate diagonalisation of $`\widehat{\rho }`$ is this coordinate basis.

It is important to note that this high-temperature approximation is well defined after a time-scale of the order of $`1/(k_\mathrm{B}T)\gamma _0/D`$ (with $`\mathrm{}=1`$). The relevant regime of evolution for our systems takes place at times comfortably larger than this time-scale.

### B Numerical Analysis

We have solved equation (21) numerically, in the high temperature limit, for different values of the diffusion coefficient $`D`$, in order to illustrate its relevance in the quantum to classical transition. Details are given elsewhere. We have chosen as initial condition a Gaussian state centered at $`x_0=p_0=0`$ with minimal uncertainty ($`\sigma _x^2=0.5`$ and $`\sigma _p^2=0.5`$). The Wigner function is initially positive definite, and different from zero only near the top of the potential. We have set the coupling constant $`\lambda =0.01`$, the renormalized frequencies $`\stackrel{~}{\mathrm{\Omega }}=\omega _n=1`$ (we are measuring time in units of $`\stackrel{~}{\mathrm{\Omega }}`$) and the bare masses also equal to one.

It is illustrative to examine first the exact result when the environment is absent (for this case we have solved numerically the Schrödinger equation). The initially Gaussian Wigner function begins to squeeze in the $`x=p`$ direction and, before the spinodal time ($`t_{\mathrm{sp}}2.3`$) it becomes a non-positive function (Fig. 1). Our definition of the spinodal time is that time at which $`\widehat{x}^2_t=2\mathrm{\Omega }_0^2/\lambda `$, its symmetry-breaking value. During the evolution, the Wigner function covers all the phase space (Fig. 2) and it is clear that it is not possible to consider it as a classical probability distribution. Although we started with a special initial state (Gaussian with minimum uncertainty), the non-linearities of the potential make the Wigner function a non-positive distribution.

Let us now consider a coupling with an environment such that the normal diffusion coefficient is $`D=0.01`$. As expected, the evolution of the Wigner function is similar to the previous one at early times (Figs. 3 and 4). However, as can be seen from Figs. 5 and 6, at long times it becomes positive definite and peaked around the classical phase space.

The effect of the environment is more dramatic for larger values of the diffusion coefficient. In our last example, $`D=1`$, the quantum to classical transition takes place almost instantaneously, even before the quantum particle pass through the spinodal point (Figs. 7 and 8).

It is interesting to note that, as the diffusion coefficient grows, the amplitude of the Wigner function decreases. This is due to the fact that the decoherence increases with $`D`$. The reduced density matrix diagonalizes. As a consequence, its ‘Fourier transform’, the reduced Wigner function, spreads out.

Our numerical results show explicitly that the existence of the environment is crucial in the quantum to classical transition. The decoherence time depends on the temperature and the coupling between system and environment through the diffusion coefficient $`D`$ .

## III Phase transitions in field theory

We can now tackle the full quantum $`\varphi `$-field theory of (1). A simplified analysis of the model has been given elsewhere. As before, the onset of classical behaviour is due to the environment. For the infinite degree of freedom quantum field, the field ordering after the transition begins is due to the growth in amplitude of its unstable long-wavelength modes. For these modes the environment consists of the short-wavelength modes of the field, together with all the other fields $`\chi _\mathrm{a}`$ with which the $`\varphi `$ inescapably interacts. The inclusion of explicit environment fields $`\chi _\mathrm{a}`$ is both a reflection of the fact that a scalar field in isolation is physically unrealistic, as well as providing us with a systematic approximation scheme.

Although the system field can never avoid the decohering environment of its short-wavelength modes, to demonstrate the effect of an environment we first consider the case in which it is taken to be composed only of the fields $`\chi _\mathrm{a}`$. The short-wavelength modes of the $`\varphi `$ field will be considered afterwards. Specifically, we take the simplest scalar classical action

$$S[\varphi ,\chi ]=S_{\mathrm{syst}}[\varphi ]+S_{\mathrm{env}}[\chi ]+S_{\mathrm{int}}[\varphi ,\chi ],$$
(22)
$`S_{\mathrm{syst}}[\varphi ]={\displaystyle d^4x\left\{\frac{1}{2}_\mu \varphi ^\mu \varphi +\frac{1}{2}\mu ^2\varphi ^2\frac{\lambda }{4!}\varphi ^4\right\}},`$ (23)
$`S_{\mathrm{env}}[\chi _\mathrm{a}]={\displaystyle \underset{\mathrm{a}=1}{\overset{N}{}}}{\displaystyle d^4x\left\{\frac{1}{2}_\mu \chi _\mathrm{a}^\mu \chi _\mathrm{a}\frac{1}{2}m_\mathrm{a}^2\chi _\mathrm{a}^2\right\}},`$ (24)
$`S_{\mathrm{int}}[\varphi ,\chi ]={\displaystyle \underset{\mathrm{a}=1}{\overset{N}{}}}{\displaystyle \frac{g_\mathrm{a}}{8}}{\displaystyle d^4x\varphi ^2(x)\chi _\mathrm{a}^2(x)},`$ (25)
where $`\mu ^2,m^2>0`$.

We see that the quantum mechanics action of (5) is, in large part, a simplified version of (22). To make any progress it is important that the environment be as simple as possible. To that end we have included no explicit self-interaction of the $`\chi _\mathrm{a}`$ fields, which interact through the $`\varphi `$ field as intermediary. However, there is one significant difference in that the interactions of the $`\varphi `$-field with the environment in (25) are quadratic, and not linear as in (5). The more conventional choice of a linear coupling to the environment was made in the previous section (and ) to give a model directly comparable to similar particle models with no symmetry breaking, for which much work has been done. Although it has been adopted for quantum field theories a linear-linear coupling is inappropriate and we have taken the simplest quadratic-quadratic couplings in $`S_{\mathrm{int}}`$ for convenience. The extension to Yukawa couplings is straightforward, and will be given elsewhere. We are reminded that, for non-linear couplings like $`x^\mathrm{n}q_\mathrm{i}^\mathrm{m}`$ in the quantum mechanical problem, one expects the master equation to contain terms of the form $`iD^{(\mathrm{n},\mathrm{m})}(t)(x^\mathrm{n}x^\mathrm{n})^2\rho _\mathrm{r}`$. We shall find a similar effect here. The non-linear coupling to the environment is crucial to our conclusions.

Further, in our present model, the environment has a strong impact upon the system-field, but not vice-versa, whenever possible. The simplest way to implement this is to take a large number $`N1`$ of $`\chi _\mathrm{a}`$ fields with comparable masses $`m_\mathrm{a}\mu `$ weakly coupled to the $`\varphi `$, with $`\lambda `$; $`g_\mathrm{a}1`$. Thus, at any step, there are $`N`$ weakly coupled environmental fields influencing the system field, but only one system field to back-react upon the explicit environment.

### A Initial conditions

Before we examine the model in detail, there are some general observations to be made about initial conditions, and the way in which the transition is implemented. Like any simple scalar theory the model displays a continuous transition at a temperature $`T_\mathrm{c}`$,

$`T_\mathrm{c}^2={\displaystyle \frac{\mu ^2}{\lambda +g_\mathrm{a}}}\mu ^2`$

The environmental fields $`\chi _\mathrm{a}`$ reduce $`T_\mathrm{c}`$ and, in order that $`T_\mathrm{c}^2\mu ^2`$, we must take

$`\lambda +{\displaystyle g_\mathrm{a}}1.`$

For order of magnitude estimations it is sufficient to take identical $`g_\mathrm{a}=\overline{g}/\sqrt{N}`$, whereby

$`11/\sqrt{N}\overline{g}\lambda .`$

i.e. the $`\chi _\mathrm{a}`$ ’tadpole’ diagrams completely overwhelm the $`\varphi `$ self-interaction tadpole diagram in generating the $`\varphi `$ thermal mass.

With $`\eta =\sqrt{6\mu ^2/\lambda }`$ determining the position of the minima of the potential and the final value of the order parameter, this choice of coupling and environments gives the hierarchy of scales

$`\mu ^2T_\mathrm{c}^2=O\left({\displaystyle \frac{\eta ^2}{\sqrt{N}}}\right)\eta ^2,`$

important in establishing a reliable approximation scheme. Further, with this choice the dominant hard loop contribution of the $`\varphi `$-field to the $`\chi _\mathrm{a}`$ thermal masses is

$`\delta m_T^2=O(\overline{g}T_\mathrm{c}^2/\sqrt{N})=O(\mu ^2/N)\mu ^2.`$

Similarly, the two-loop (setting sun) diagram which is the first to contribute to the discontinuity of the $`\chi `$-field propagator is of magnitude

$`\overline{g}^2T_\mathrm{c}^2/N=O(g\mu ^2/N^{3/2})\delta m_T^2,`$

in turn. That is, the effect of the thermal bath on the propagation of the environmental $`\chi `$-fields is ignorable.

This was our intention in model-making; to construct an environment that reacted on the system field, but was not reacted upon by it to any significant extent. We stress that this is not a Hartree or large-N approximation of the type that, to date, has been the standard way to proceed for a closed system.

We shall assume that the initial states of the system and environment are both thermal, at a temperature $`T_0>T_\mathrm{c}`$. We then imagine a change in the global environment (e.g. expansion in the early universe) that can be characterised by a change in temperature from $`T_0`$ to $`T_\mathrm{f}<T_\mathrm{c}`$. That is, we do not attribute the transition to the effects of the environment-fields. On incorporating the hard thermal loop tadpole diagrams of the $`\chi `$ (and $`\varphi `$) fields in the $`\varphi `$ mass term leads to the effective action for $`\varphi `$ quasiparticles,

$`S_{\mathrm{syst}}^{\mathrm{eff}}[\varphi ]={\displaystyle d^4x\left\{\frac{1}{2}_\mu \varphi ^\mu \varphi \frac{1}{2}m_\varphi ^2(T_0)\varphi ^2\frac{\lambda }{4!}\varphi ^4\right\}}`$

where $`m_\varphi ^2(T_0)=\mu ^2(1T_0^2/T_c^2)=M^2>0`$. As a result, we can take an initial factorised density matrix at temperature $`T_0`$ of the form $`\widehat{\rho }[T_0]=\widehat{\rho }_\varphi [T_0]\widehat{\rho }_\chi [T_0]`$, where $`\widehat{\rho }_\varphi [T_0]`$ is determined by the quadratic part of $`S_{\mathrm{syst}}^{\mathrm{eff}}[\varphi ]`$ and $`\widehat{\rho }_\chi [T_0]`$ by $`S_{\mathrm{env}}[\chi _\mathrm{a}]`$. Yet again, the many $`\chi _\mathrm{a}`$ fields have a large effect on $`\varphi `$, but the $`\varphi `$-field has negligible effect on the $`\chi _\mathrm{a}`$. Provided the change in temperature is not too slow the exponential instabilities of the $`\varphi `$-field grow so fast that the field has populated the degenerate vacua well before the temperature has dropped to zero. Since the temperature $`T_\mathrm{c}`$ has no particular significance for the environment field, for these early times we can keep the temperature of the environment fixed at $`T_\chi =T_0=𝒪(T_\mathrm{c})`$ (our calculations are only at the level of orders of magnitude).

Since it is the system-field $`\varphi `$ field whose behaviour changes dramatically on taking $`T_\varphi `$ through $`T_\mathrm{c}`$, in this paper we adopt an instantaneous quench for $`T_\varphi `$ from $`T_0`$ to $`T_\mathrm{f}=0`$ at time $`t=0`$, in which $`m_\varphi ^2(T)`$ changes sign and magnitude instantly, concluding with the value $`m_\varphi ^2(t)=\mu ^2,t>0`$ that we cited in (1). An instantaneous quench is sufficient to demonstrate the rapidity with which the environment forces the system field to become classical. Meanwhile, for simplicity the $`\chi _\mathrm{a}`$ masses are fixed at the common value $`m\mu `$.

### B Tracing out the $`\chi `$ fields

At time $`t>0`$ the reduced density matrix $`\rho _\mathrm{r}[\varphi ^+,\varphi ^{},t]=\varphi ^+|\widehat{\rho }_r(t)|\varphi ^{}`$ is now

$`\rho _\mathrm{r}[\varphi ^+,\varphi ^{},t]={\displaystyle 𝒟\chi _\mathrm{a}\rho [\varphi ^+,\chi _\mathrm{a},\varphi ^{},\chi _\mathrm{a},t]},`$

where $`\rho [\varphi ^+,\chi _\mathrm{a}^+,\varphi ^{},\chi _\mathrm{a}^{},t]=\varphi ^+\chi _\mathrm{a}^+|\widehat{\rho }(t)|\varphi ^{}\chi _\mathrm{a}^{}`$ is the full density matrix. Since we would like to be able to distinguish between different classical system-field configurations evolving after the transition, we will only interested in the field-configuration basis for this reduced density matrix (in analogy with the quantum Brownian motion model of the previous section). The environment will have had the effect of making the system essentially classical once $`\rho _\mathrm{r}(t)`$ is, effectively, diagonal. \[Our earlier comments on dephasing remain valid, and we stress again that our understanding of what constitutes classical behaviour is essentially different from the dephasing effects found in Refs. ).\] Quantum interference can then be ignored and the system is said to have decohered. At the same time, we obtain a classical probability distribution from the diagonal part of $`\rho _\mathrm{r}(t)`$, or equivalently, by means of the reduced Wigner functional. For weak coupling there will be no ’recoherence’ at later times in which the sense of classical probability will be lost.

Its temporal evolution is given by

$`\rho _\mathrm{r}[\varphi _\mathrm{f}^+,\varphi _\mathrm{f}^{},t]={\displaystyle 𝑑\varphi _\mathrm{i}^+𝑑\varphi _\mathrm{i}^{}J_\mathrm{r}[\varphi _\mathrm{f}^+,\varphi _\mathrm{f}^{},t|\varphi _\mathrm{i}^+,\varphi _\mathrm{i}^{},t_0]\rho _\mathrm{r}[\varphi _\mathrm{i}^+\varphi _\mathrm{i}^{},t_0]},`$

where $`J_\mathrm{r}`$ is the reduced evolution operator

$$J_\mathrm{r}[\varphi _\mathrm{f}^+,\varphi _\mathrm{f}^{},t|\varphi _\mathrm{i}^+,\varphi _\mathrm{i}^{},t_0]=_{\varphi _\mathrm{i}^+}^{\varphi _\mathrm{f}^+}𝒟\varphi ^+_{\varphi _\mathrm{i}^{}}^{\varphi _\mathrm{f}^{}}𝒟\varphi ^{}e^{i\{S[\varphi ^+]S[\varphi ^{}]\}}F[\varphi ^+,\varphi ^{}].$$
(26)
The Feynman-Vernon  influence functional $`F[\varphi ^+,\varphi ^{}]`$ is defined as

$`F[\varphi ^+,\varphi ^{}]={\displaystyle 𝑑\chi _{\mathrm{ai}}^+𝑑\chi _{\mathrm{ai}}^{}\rho _\chi [\chi _{\mathrm{ai}}^+,\chi _{\mathrm{ai}}^{},t_0]𝑑\chi _{\mathrm{af}}}`$ (27)
$`\times `$ $`{\displaystyle _{\chi _{\mathrm{ai}}^+}^{\chi _{\mathrm{af}}}}𝒟\chi _\mathrm{a}^+{\displaystyle _{\chi _{\mathrm{ai}}^{}}^{\chi _{\mathrm{af}}}}𝒟\chi _\mathrm{a}^{}\mathrm{exp}(i\{S_{\mathrm{env}}[\chi _\mathrm{a}^+]+S_{\mathrm{int}}[\varphi ^+,\chi _\mathrm{a}^+])`$ (28)
$`\times `$ $`\mathrm{exp}\left(i\{S_{\mathrm{env}}[\chi _\mathrm{a}^{}]+S_{\mathrm{int}}[\varphi ^{},\chi _\mathrm{a}^{}]\}\right).`$ (29)
Beginning from this initial distribution, peaked around $`\varphi =0`$, we follow the evolution of the system under the influence of the environment fields, with Hamiltonian determined from (22). From the influence functional we define the influence action $`\delta A[\varphi ^+,\varphi ^{}]`$ as

$$F[\varphi ^+,\varphi ^{}]=\mathrm{exp}i\delta A[\varphi ^+,\varphi ^{}].$$
(30)
The total coarse-grained effective action is

$`A[\varphi ^+,\varphi ^{}]=S[\varphi ^+]S[\varphi ^{}]+\delta A[\varphi ^+,\varphi ^{}].`$

We will calculate the influence action to one loop (two vertices) for $`N`$ large using closed-time path correlators. It is the imaginary part which contains the information about the onset of classical behaviour. If $`\mathrm{\Delta }=\frac{1}{2}(\varphi ^{+2}\varphi ^2)`$, we find

$$\mathrm{Im}\delta A=\frac{\overline{g}^2}{16}d^4xd^4y\mathrm{\Delta }(x)N_\mathrm{q}(x,y)\mathrm{\Delta }(y).$$
(31)
In (31) $`N_\mathrm{q}(xy)=\mathrm{Re}G_{++}^2(x,y)`$ is the noise (diffusion) kernel, where $`G_{++}`$ is the relevant closed-timepath correlator of the $`\chi `$-field at temperature $`T_0`$. Non-leading one-loop terms are smaller by a factor $`𝒪(N^{1/2})`$.

The first step in the evaluation of the master equation is the calculation of the density matrix propagator $`J_\mathrm{r}`$ from Eq. (26). In order to estimate the functional integration which defines the reduced propagator, we perform a saddle point approximation

$`J_\mathrm{r}[\varphi _\mathrm{f}^+,\varphi _\mathrm{f}^{},t|\varphi _\mathrm{i}^+,\varphi _\mathrm{i}^{},t_0]\mathrm{exp}iA[\varphi _{\mathrm{cl}}^+,\varphi _{\mathrm{cl}}^{}],`$

where $`\varphi _{\mathrm{cl}}^\pm `$ is the solution of the equation of motion $`\frac{\delta ReA}{\delta \varphi ^+}|_{\varphi ^+=\varphi ^{}}=0`$ with boundary conditions $`\varphi _{\mathrm{cl}}^\pm (t_0)=\varphi _\mathrm{i}^\pm `$ and $`\varphi _{\mathrm{cl}}^\pm (t)=\varphi _\mathrm{f}^\pm `$. It is very difficult to solve this equation analytically. For simplicity, we assume that the system-field contains only one Fourier mode with $`\stackrel{}{k}=\stackrel{}{k}_0`$. We are motivated in this by the observation that the long-wavelength modes start growing exponentially as soon as the quench is performed and rapidly dominate the fluctuation power spectrum. Modes with $`|k_0|^2>\mu ^2`$ will oscillate.

Further, we are interested primarily in the order-parameter

$`\varphi (t)=\varphi (k_0=0,t)=\underset{V\mathrm{}}{lim}{\displaystyle \frac{1}{V}}{\displaystyle _{𝐱V}}d^3x\varphi (𝐱,t),`$

and all our subsequent calculations will be for theis zero-frequency mode. We write the spatially-constant classical solution as $`\varphi (s)=f(s,t)`$ where $`f(s,t)`$ satisfies the boundary conditions $`f(0,t)=\varphi _\mathrm{i}`$ and $`f(t,t)=\varphi _\mathrm{f}`$. Qualitatively, $`f(s,t)`$ grows exponentially with $`s`$ for $`tt_{\mathrm{sp}}`$, and oscillates for $`t_{\mathrm{sp}}<s<t`$ when $`t>t_{\mathrm{sp}}`$. We shall therefore approximate its time dependence for $`tt_{\mathrm{sp}}`$ as

$$f(s,t)=\varphi _\mathrm{i}u_1(s,t)+\varphi _\mathrm{f}u_2(s,t),$$
(32)
where $`u_1(0,t)=1`$, $`u_1(t,t)=0`$ and $`u_2(0,t)=0`$, $`u_2(t,t)=1`$, with solution

$`u_1(s,t)={\displaystyle \frac{\mathrm{sinh}[\mu (ts)]}{\mathrm{sinh}(\mu t)}},u_2(s,t)={\displaystyle \frac{\mathrm{sinh}(\mu s)}{\mathrm{sinh}(\mu t)}}.`$

We shall justify the use of the linear equation later. In fact, we shall see then that it is not an unreasonable approximation until almost $`t_{\mathrm{sp}}`$.

In order to solve the master equation we must compute the final time derivative of the propagator $`J_\mathrm{r}`$, and after that eliminate the dependence on the initial field configurations $`\varphi _\mathrm{i}^\pm `$ coming from the classical solutions $`\varphi _{\mathrm{cl}}^\pm `$. This is the usual procedure, see Ref..

As we are solely interested in the onset of classical behaviour, it is sufficient to calculate the correction to the usual unitary evolution coming from the noise kernel. For clarity we drop the suffix $`\mathrm{f}`$ on the final state fields. If $`\mathrm{\Delta }=(\varphi ^{+2}\varphi ^2)/2`$ for the final field configurations, then the master equation for $`\rho _r(\varphi ^+,\varphi ^{},t)`$ is

$$i\dot{\rho }_\mathrm{r}=\varphi ^+|[H,\widehat{\rho }_\mathrm{r}]|\varphi ^{}i\frac{\overline{g}^2}{16}V\mathrm{\Delta }^2D(t)\rho _\mathrm{r}+\mathrm{}$$
(33)
The presence of $`\mathrm{\Delta }`$ in (33), rather than $`\varphi _+\varphi _{}`$, is a consequence of the quadratic coupling to the environment in $`S_{\mathrm{int}}`$. Since our main concern is with the diagonalisation of $`\rho _\mathrm{r}`$, it is not necessary to consider the evolution equation for the Wigner functional, as has been shown in the literature of quantum Brownian motion . The volume factor $`V`$ that appears in the master equation is due to the fact we are considering a density matrix which is a functional of two different field configurations, $`\varphi ^\pm (\stackrel{}{x})=\varphi ^\pm `$, which are spread over all space. The time dependent diffusion coefficient $`D_\chi (t)`$ due to each of the many external environmental $`\chi `$ fields is given by

$$D_\chi (t)=3_0^t𝑑su(s,t)\mathrm{Re}G_{++}^2(0;ts),$$
(34)
where

$$u(s,t)=\left[u_2(s,t)\frac{\dot{u}_2(t,t)}{\dot{u}_1(t,t)}u_1(s,t)\right]^2.$$
(35)
For the case in hand of an instantaneous quench, $`u(s,t)=\mathrm{cosh}^2\mu (ts)`$ when $`tt_{\mathrm{sp}}`$, and is an oscillatory function of time when $`t>t_{\mathrm{sp}}`$.

Although $`G_{++}`$ is oscillatory at all times, the exponential growth of $`u(t)`$ enforces a similar behaviour on $`D_\chi (t)`$. In the high temperature limit ($`k_\mathrm{B}T\mu `$), the explicit expression for that contribution to the diffusion coefficient due to the $`\chi _\mathrm{a}`$ fields alone is

$$D_\chi (t)=\frac{5(k_\mathrm{B}T)^2}{64\pi ^2}_0^t𝑑su(s,t)_0^{\mathrm{}}𝑑p\frac{p^2}{\omega ^4}\mathrm{cos}[2\omega s],$$
(36)
where $`\omega ^2=p^2+m^2`$ and we have set $`m^2=\mu ^2`$.

For times $`\mu t1`$, the integration in (36) is dominated by the behaviour at $`s=0`$:

$$D_\chi (t)(k_\mathrm{B}T_0/\mu )^2u(0,t)(k_\mathrm{B}T_0/\mu )^2\mathrm{exp}[2\mu t].$$
(37)
The spinodal time $`t_{\mathrm{sp}}`$ is again defined as the time for which $`\varphi ^2_t\eta ^2=6\mu ^2/\lambda `$. For $`t>t_{\mathrm{sp}}`$ the diffusion coefficient stops growing, and oscillates around $`D(t=t_{\mathrm{sp}})`$.

### C Short wavelength modes

In our present model the environment fields $`\chi _\mathrm{a}`$ are not the only decohering agents. The environment is also constituted by the short-wavelength modes of the self-interacting field $`\varphi `$. Therefore, we split the field as $`\varphi =\varphi _<+\varphi _>`$, where the system-field $`\varphi _<`$ contains the modes with wavelengths longer than the critical value $`\mu ^1`$, while the bath or environment-field contains wavelengths shorter than $`\mu ^1`$. This gives an additional one-loop contribution $`D_\varphi (t)`$ to the diffusion function with the same $`u(s,t)`$ but a $`G_{++}`$ constructed form the short-wavelength modes of the $`\varphi `$-field as it evolves from the top of the potential hill. Without the additional powers of $`N^1`$ to order contributions summation of loop diagrams is essential to get a reliable $`G_{++}`$. However, it is not necessary to calculate $`D_\varphi (t)`$ in order to get a good estimate of $`t_D`$. Since the contribution of $`D_\varphi (t)`$ to the overall diffusion function is positive we can derive an upper bound on the decoherence time $`t_D`$ from the reliable diffusion functions $`D_\chi `$ alone.

However, we would not expect the inclusion of the $`\varphi `$-field to give a qualitative change. Specifically, we note that, if we take $`\varphi `$-field propagators dressed by only the simplest tadpole diagrams, and ignore two-loop self-interaction diagrams then the diffusion correction due to the $`\varphi `$ loop is now similar to that of the $`\chi `$ loops. The effect is that the short-wavelength modes in the one-loop diagrams from which they are calculated have been kept at the initial temperature $`T_0`$, on the grounds that passing through the transition quickly has no effect on them. The quench mimics a slower evolution of temperature in which only the long-wavelength modes show instabilities that the transition induces. That is, with $`\overline{g}\lambda `$ and no $`1/N`$ factor, the short-wavelength $`\varphi `$ modes would have the same effect on the dissipation, qualitatively, as all the explicit environmental fields put together. However, at an order of magnitude level there is no change, since the effect is to replace $`\overline{g}^2`$ by $`\overline{g}^2+O(\lambda ^2)=O(\overline{g}^2)`$.

### D The decoherence time

Using the positivity of $`D_\varphi `$ we estimate the decoherence time $`t_\mathrm{D}`$ for the model by considering the approximate solution to the master equation (33),

$`\rho _\mathrm{r}[\varphi _<^+,\varphi _<^{};t]\rho _\mathrm{r}^\mathrm{u}[\varphi _<^+,\varphi _<^{};t]\mathrm{exp}\left[V\mathrm{\Gamma }{\displaystyle _0^t}𝑑sD_\chi (s)\right],`$

where $`\rho _\mathrm{r}^\mathrm{u}`$ is the solution of the unitary part of the master equation (i.e. without environment). It is obvious from this (and also from (33)), that the diagonal density matrix just evolves like the unitary matrix (the environment has almost no effect on the diagonal part of $`\rho _\mathrm{r}`$). In terms of the dimensionless fields $`\overline{\varphi }=(\varphi _<^++\varphi _<^{})/2\mu ,`$ and $`\delta =(\varphi _<^+\varphi _<^{})/2\mu `$, we have $`\mathrm{\Gamma }=(1/16)\overline{g}^2\mu ^4\overline{\varphi }^2\delta ^2`$.

The system behaves classically when $`\rho _\mathrm{r}`$ is appropriately diagonal. We therefore look at the ratio

$$\left|\frac{\rho _\mathrm{r}[\overline{\varphi }+\delta ,\overline{\varphi }\delta ;t]}{\rho _\mathrm{r}[\overline{\varphi },\overline{\varphi };t]}\right|\left|\frac{\rho _\mathrm{r}^\mathrm{u}[\overline{\varphi }+\delta ,\overline{\varphi }\delta ;t]}{\rho _\mathrm{r}^\mathrm{u}[\overline{\varphi },\overline{\varphi };t]}\right|\mathrm{exp}\left[V\mathrm{\Gamma }_0^t𝑑sD_\chi (s)\right].$$
(38)
It is not possible to obtain an analytic expression for the ratio of unitary density matrices that appears in Eq.(38). The simplest approximation is to neglect the self-coupling of the system field . In this case the unitary density matrix remains Gaussian at all times as

$$\left|\frac{\rho _\mathrm{r}^\mathrm{u}[\overline{\varphi }+\delta ,\overline{\varphi }\delta ;t]}{\rho _\mathrm{r}^\mathrm{u}[\overline{\varphi },\overline{\varphi };t]}\right|=\mathrm{exp}[\frac{T_\mathrm{c}}{\mu }\delta ^2p^1(t)],$$
(39)
where $`p^1(t)`$, essentially $`\mu ^2\varphi ^2_t^1`$, decreases exponentially with time to a value $`𝒪(\lambda )`$. This approximation can be improved by means of a Hartree-like approximation . In this case the ratio is still given by Eq.(39), but now $`p^1(t)`$ decreases more slowly as $`t`$ approaches $`t_{\mathrm{sp}}`$. In any case, in the unitary part of the reduced density matrix the non-diagonal terms are not suppressed. Therefore, in order to obtain classical behaviour, the relevant part of the reduced density matrix is the term proportional to the diffusion coefficient in Eq.(38), since it is this that enforces its diagonalisation.

The decoherence time $`t_\mathrm{D}`$ sets the scale after which we have a classical system-field configuration. According to our previous discussion, it can be defined as the solution to

$$1V\mathrm{\Gamma }_0^{t_D}𝑑sD(s).$$
(40)
It corresponds to the time after which we are able to distinguish between two different field amplitudes (inside a given volume $`V`$).

Suppose we reduce the couplings $`\overline{g}`$, $`\lambda `$ of the system $`\varphi `$-field to its environment. Since, as a one-loop construct, $`\mathrm{\Gamma }\overline{g}^2,\lambda ^2`$ our first guess would be that as $`\overline{g}`$, $`\lambda `$ decrease, then $`t_D`$ increases and the system takes longer to become classical. This is not really the case. The reason is twofold. Firstly, there is the effect that $`\mathrm{\Gamma }T_0^2`$, and $`T_0^2\lambda ^1`$ is non-perturbatively large for a phase transition. Secondly, because of the non-linear coupling to the environment, obligatory for quantum field theory, $`\mathrm{\Gamma }\overline{\varphi }^2`$. The completion of the transition finds $`\overline{\varphi }^2\eta ^2\lambda ^1`$ also non-perturbatively large. This suggests that $`\mathrm{\Gamma }`$, and hence $`t_D`$, can be independent of $`\lambda `$. The situation would be different for a linear coupling to the environment, for which $`\overline{\varphi }^2`$ would not be present, or a cold initial state in which $`\varphi `$ is peaked about $`\varphi =0`$. In fact, the situation is a little more complicated, but the corollary that $`t_D`$ does not exceed $`t_{\mathrm{sp}}`$ as the couplings become weaker remains true.

In order to quantify the decoherence time we have to fix the values of $`V`$, $`\delta `$, and $`\overline{\varphi }`$. $`V`$ is understood as the minimal volume inside which we do not accept coherent superpositions of macroscopically distinguishable states for the field. Thus, our choice is that this volume factor is $`𝒪(\mu ^3)`$ since $`\mu ^1`$ (the Compton wavelength) is the smallest scale at which we need to look. Inside this volume, we do not discriminate between field amplitudes which differ by $`𝒪(\mu )`$, and therefore we take $`\delta 𝒪(1)`$. For $`\overline{\varphi }`$ we set $`\overline{\varphi }^2𝒪(\alpha /\lambda )`$, where $`\lambda \alpha 1`$ is to be determined self-consistently.

Note that the diagonalisation of $`\rho _t`$ occurs in time as an exponential of an exponential. As a result, decoherence occurs extremely quickly, but not so quickly that $`\mu t1`$. Consequently, in order to evaluate the decoherence time in our model, we have to use Eq.(37). On taking the equality in (40) we find

$$\mathrm{exp}[2\mu t_D]\frac{\lambda \sqrt{N}\overline{g}}{\overline{g}^2\alpha }=𝒪(\frac{\sqrt{N}}{\alpha }),$$
(41)
whereby

$`\mu t_D{\displaystyle \frac{1}{4}}\mathrm{ln}N{\displaystyle \frac{1}{2}}\mathrm{ln}\alpha \mathrm{ln}({\displaystyle \frac{\eta }{T_\mathrm{c}\sqrt{\alpha }}}).`$

This is an upper bound on $`t_D`$, but probably still qualitatively correct. The value of $`\alpha `$ is determined as $`\alpha \sqrt{\mu /T_\mathrm{c}}`$ from the condition that, at time $`t_D`$, $`|\varphi |^2_t\alpha \eta ^2`$. Since $`\alpha 1`$, in principle, the field has not diffused far from the top of the hill before it is behaving classically.

For comparison, we find $`t_{\mathrm{sp}}`$, for which $`\varphi ^2_t\eta ^2`$, given by

$$\mathrm{exp}[2\mu t_{\mathrm{sp}}]𝒪(\frac{\eta ^2}{\mu T_\mathrm{c}}).$$
(42)
The exponential factor, as always, arises from the growth of the unstable long-wavelength modes. The factor $`T_\mathrm{c}^1`$ comes from the $`\mathrm{coth}(\beta \omega /2)`$ factor that encodes the initial Boltzmann distribution at temperature $`TT_\mathrm{c}`$. As a result,

$$\mu t_{\mathrm{sp}}\mathrm{ln}(\frac{\eta }{\sqrt{\mu T_\mathrm{c}}}),$$
(43)
whereby $`1<\mu t_D\mu t_{\mathrm{sp}}`$, with

$$\mu t_{\mathrm{sp}}\mu t_D\frac{1}{4}\mathrm{ln}(\frac{T_\mathrm{c}}{\mu })>1,$$
(44)
for weak enough coupling, or high enough initial temperatures (we have taken $`T_0T_\mathrm{c}`$ throughout). This is our main result, that for the physically relevant modes (with small $`k_0`$) classical behaviour has been established before the spinodal time, when the ground states have became populated.

### E Back-reaction

We can now justify our earlier assumption that, for an instantaneous quench, nonlinear behaviour only becomes important just before the spinodal time. To see this, we adopt the Hartree approximation, in which the equations of motion are linearised so that the field modes $`f_k^\pm `$ now satisfy the equation

$`\left[{\displaystyle \frac{d^2}{dt^2}}+𝐤^2\mu ^2(t)\right]f_k^\pm (t)=0,`$

where

$$\mu ^2(t)=m^2(t)+4\lambda \frac{d^3p}{(2\pi )^3}C(p)[f_p^+(t)f_p^{}(t)1]$$
(45)
and $`C(p)=(\mathrm{coth}\beta \omega /2)/2\omega `$, $`\omega ^2=p^2+m^2`$.

A similar result would be obtained by extending our initial $`O(2)`$ theory to an $`O(M)`$ theory in the large-M limit before taking $`M=2`$.

Our coarse-graining retains only the unstable modes in the integral, which suggests the approximate hybrid self-consistent equation for $`\mu ^2(t)`$,

$$\mu ^2(t)\mu ^2C\lambda \frac{T\mu }{(\mu t_{\mathrm{sp}})^{3/2}}\mathrm{exp}\left(2_0^t𝑑t^{}\mu (t^{})\right),$$
(46)
($`C=O(1)`$), which has the exponential growth of the WKB solution, but non-singular behaviour when $`\mu (t)0`$. The exact solution to Eq.(46) for $`t<t_{\mathrm{sp}}`$ is $`\mu (t)=\mu \mathrm{tanh}\mu (t_{\mathrm{sp}}t)`$, irrespective of the values of the temperature $`T`$ and the coupling strength. In fact, we anticipated this in (43), when we estimated $`t_{\mathrm{sp}}`$ on the assumption that the backreaction would only take effect in the final moments.

That is, the theory only ceases to behave like a free Gaussian theory with upside-down potential at a time $`t_\mathrm{B}`$,

$$t_{\mathrm{sp}}t_\mathrm{B}=O(\mu ^1).$$
(47)
It follows that $`t_\mathrm{B}t_D`$ in our ordering of scales $`T_\mathrm{c}\mu `$ but, in practice $`T_\mathrm{c}`$ needs to be at least an order of magnitude larger than $`\mu `$ for this to be the case.

When (44) is valid, we see that $`\rho _\mathrm{r}`$ becomes diagonal before non-linear terms could be relevant. In this sense, classical behaviour has been achieved before quantum effects could destroy the positivity of the Wigner function $`W_\mathrm{r}`$.

## IV Final remarks

We have shown that, in our model, $`\rho _\mathrm{r}`$ becomes diagonal before the spinodal time at which the order parameter field has first populated the ground-state values of the theory. Further, it can also be diagonal before non-linear terms are relevant. That is, decoherence can be achieved before quantum effects destroy the positivity of the Wigner function $`W_\mathrm{r}`$. Really, our $`t_\mathrm{D}`$ sets the time after which we have a classical probability distribution (positive definite) even for times $`t>t_{\mathrm{sp}}`$. The existence of the environment is crucial in doing this. Of course, for non-Gaussian or delocalised (in the field space) initial states, it is clear that $`W_\mathrm{r}`$ will be non-positive definite even in the linear regime, and therefore $`t_\mathrm{D}`$ should be smaller than the one we evaluated here. In the present work, $`t_\mathrm{D}`$ is the classicalisation time-bound.

This result goes in the direction of justifying the use of classical numerical simulations for the analysis of the dynamics of the long-wavelengths modes after the quench.

###### Acknowledgements.

F.C.L. and F.D.M. were supported by Universidad de Buenos Aires, CONICET (Argentina), Fundación Antorchas, and ANPCyT. R.J.R. was supported in part by the COSLAB programme of the European Science Foundation. We also thank the organisers of the Peyresq VI meeting.


text:
# Time Domain Computation of a Nonlinear Nonlocal Cochlear Model with Applications to Multitone Interaction in Hearing

## 1 Introduction

Voice signal processing algorithms have received increasing attention in recent years for improving the design of hearing devices and for evaluating acoustic theories of peripheral auditory systems, see Meddis et al  among others. A fundamental issue a computational method needs to resolve is the nonlinear interaction of acoustic waves of different frequencies and the resulting tonal suppression effects. A nonlinear filter bank approach is developed in , based on knowledge of auditory responses. Nonlinearities are known to originate in the cochlea and are further modified in higher level auditory pathways. The cochlear mechanics has first principle descriptions, and so partial differential equations (PDEs) form a natural mathematical framework to initiate computation. However, in vivo cochlear dynamics is not a pure mechanical problem, and neural couplings are present to modify responses. To incorporate both of these aspects, we aim to develop a first-principle based PDE approach to voice signal processing, where the neural aspect is introduced in the model phenomenologically.

Cochlear modeling has had a long history, and various linear models have been studied at length by analytical and numerical methods, see Keller and Neu , Leveque, Peskin and Lax , Lighthill , and references therein. It has been realized that nonlinearity is essential for multitone interactions, see  among others. The nonlinearity could be incorporated through micro-mechanics of cochlea, such as coupling of basilar membrane (BM) to inner hair cells . Or nonlinearity could be introduced phenomenologically based on spreading of electrical and neural activities between hair cells at different BM locations suggested by experimental data, see Jau and Geisler , Deng . The latter treatment turned out to be quite efficient and will be adopted in this paper. Their model with nonlinear and nonlocal BM damping will be our starting point.

Multitone interaction requires one to perform numerical simulation and analysis in the time domain as the resulting time dependent problem is strictly speaking irreducible to a steady state problem. The commonly used cochlear models, including those of  and , are intrinsically dispersive. In particular, long waves tend to propagate with little decay from entrance point (stapes) to the exit point (helicotrama). Such an issue can be avoided in linear models because a reduction to the steady state (or a frequency domain calculation using time harmonic solutions) is available by factoring out the time dependence of solutions. The dispersive phenomenon prompted us to incorporate a selective fluid longitudinal damping term in the model of Jau and Geisler , Deng , operating near the exit point. The role of such a term is to suck out the long waves accumulating near the exit. Selective positive or negative damping has been a novel way to filter images in PDE method of image processing, see Osher and Rudin  for the first work in this direction. Our treatment here is similar in spirit though on a different type of PDE.

We then present a second order semi-implicit finite difference method of our model. The method relaxes the stability time step constraint of explicit methods without going into the complexity of fully implicit schemes. As our model is nonlocal and nonlinear, semi-implicit scheme is a reasonable option for computation. Numerical examples showed that the method is efficient, and that the selective damping term indeed removes the long waves at the exit boundary point for a time domain calculation of a single tone. The examples also demonstrated that the analysis of dispersive waves and slow decay of long waves is robust, and the phenomenon persists in the nonlinear nonlocal regime of the model.

Subsequently, we show computational results on isodisplacement curves for five characteristic frequencies consistent with those of Neely  for BM displacement of 1 nm. In case of two tone interactions, our model gives both low and high side suppression, in qualitative agreement with similar studies of Geisler . For understanding three tone interactions, we fix the frequencies and amplitudes of two tones, and vary the third tone both in frequency and amplitude in a manner similar to auditory neural experiment of Deng and Geisler . Again, our model showed suppressing effect on the two tones by the third tone, in qualitative agreement with experimental data. To our best knowledge, the current paper is the first on time domain simulation of interactions of three tones using first principle based PDEs. The stage is set for further computation, analysis, and model development on more complicated tones.

The rest of the paper is organized as follows. In section 2, we present our model and its background, also analysis of dispersive properties of solutions. In section 3, we outline our semi-implicit discretization of the model. In section 4, we give model parameters, and show numerical results on dispersion of model solutions, isodisplacement curves, two and three tone interactions in comparison with existing data. The conclusions are in section 5.

## 2 Nonlinear Nonlocal Cochlear Model

### 2.1 Background and Model Setup

The cochlear modeling has had a long history driven by advancement of cochlear experiments, see  among others. A brief derivation of the one dimensional (1-D) cochlear model of transmission line (long wave) type, based on fluid mechanics and elasticity, can be found in Sondhi , see also references therein for earlier contributions, and de Boer  for later development and higher dimensional models. A general form of the 1-D model can be expressed as:

$`p_{xx}Nu_{tt}=ϵ_su_t,x(0,L),`$ (2.1)
$`p=mu_{tt}+r(x,u)u_t+s(x)u,`$ (2.2)
where $`p`$ is the fluid pressure difference across the basilar membrane (BM), $`u`$ the BM displacement, $`L`$ the longitudinal length of BM; stapes is at $`x=0`$, and helicotrema is at $`x=L`$; $`N`$ a constant depending on fluid density and cochlear channel size; $`ϵ_s`$ depends on the damping of longitudinal fluid motion ; $`m`$, $`r`$, $`s`$ are the mass, damping, and stiffness of BM per unit area, with $`m`$ a constant, $`s`$ a known function of $`x`$, and $`r`$ a function(al) of $`x`$, $`u`$.

Choosing a functional form for $`r`$ is a key step in model design. Many choices in the literature assume that $`r`$ is a local nonlinear function, . The nonlinearity is necessary to model two tone nonlinear interaction (e.g. two tone suppression effects). However, as pointed out by Kim , local form of $`r`$ is not sufficient to account for two tone suppression when suppressors are below the excitors in frequency (low side), unless a ”second filter” is incorporated, see e.g . To achieve the second filter effect on more physical grounds, micromechanical approaches have been developed to derive damping $`r`$ from motion of outer hair cells and other physiologically relevant parts of cochlear,  among others. Though this is a reasonable approach, one faces the daunting complexity of the cochlear micromechanics, see  for recent progress along this line. The alternative approach, first proposed in this context by Jau and Geisler , is phenomenological in nature. It hypothesized that there is longitudinal coupling in the cochlear partition itself, and damping $`r`$ is a convolution integral of $`u`$ with an exponential weight function on BM. Similar idea appeared also in Lyon  in coupled automatic gain control framework. Possible physiological mechanisms underlying the spatially coupled nonlocal damping were identified in Deng , and the resulting model was then used to process acoustic signals.

We shall adopt the latter nonlocal damping approach for choosing the function $`r`$. The BM damping coefficient is specified as in Jau and Geisler , Deng :

$$r(x,t,u)=r_0(x)+\gamma _0^L|u(x^{},t)|\mathrm{exp}\{|xx^{}|/\lambda \}𝑑x^{},$$
(2.3)
where $`r_0(x)`$ is the passive BM damping; $`\gamma `$ and $`\lambda `$ are positive constants. The size of $`r`$ controls how sharp the BM traveling wave envelope would be.

We shall take $`ϵ_s`$ as a function of $`x`$ so that $`ϵ_s=ϵ_s(x)0`$ is a smooth function and is supported near the vicinity of $`x=L`$, for instance:

$$ϵ_s(x)=\frac{ϵ}{1+\mathrm{exp}\{\beta (x_sx)\}},$$
(2.4)
where $`x_s`$ is a constant slightly smaller than $`L`$; $`ϵ>1`$, $`\beta >1`$, are two constants, with $`\beta `$ suitably large. Taking $`ϵ_s`$ as a spatially dependent function can be considered as a way to selectively introduce damping so that possible low frequency waves near $`x=L`$ are damped out and there is minimal wave accumulation (or reflection) close to the helicotrema (at $`x=L`$). This turns out to be an essential stabilizing mechanism for our time domain numerical computation. We will come back to this point later.

The coefficients $`m`$, $`n`$, and $`s(x)`$ are standard. The function $`s(x)`$ is based on the data of Liberman :

$$s(x)=4\pi ^2m(0.456\mathrm{exp}(\mathrm{\hspace{0.17em}4.83}(\mathrm{\hspace{0.17em}1}x/3.5))0.45)^2.$$

The physical boundary and initial conditions are:

$`p_x(0,t)=T_Mp_T(t),p(L,t)=0,`$ (2.5)
$`u(x,0)=u_t(x,0)=0,`$ (2.6)
where $`p_T(t)`$ is the input sound pressure at the eardrum; and $`T_M`$ is a bounded linear operator on the space of bounded continuous functions, with output depending on the frequency content of $`p_T(t)`$. If $`p_T=_{j=1}^{J_M}A_j\mathrm{exp}\{i\omega _jt\}+c.c.`$, a multitone input, c.c denoting complex conjugate, $`J_M`$ a positive integer, then $`T_Mp_T(t)=_{j=1}^{J_M}B_j\mathrm{exp}\{i\omega _jt\}+c.c`$, where $`B_j=a_M(\omega _j)A_j`$, c.c for complex conjugate, $`a_M()`$ a scaling function built from the filtering characteristics of the middle ear. Established data are available in Guinan and Peake . The model setup is now complete.

### 2.2 Model Properties and Analysis

We briefly go over some special regimes and solutions of the model system to illustrate the basic mathematical properties and issues. If $`m=0`$, $`r=0`$, $`ϵ=0`$, (2.1)-(2.2) reduces to a second order wave equation: $`nu_{tt}(s(x)u)_{xx}=0`$. When $`ϵ`$ is turned on, and $`\beta =0`$, we have a damped second order wave equation driven by boundary condition at $`x=0`$. Now suppose that all coefficients are constants, then (2.1)-(2.2) considered on the entire line admits planar wave solutions of the form:

$$u=u_0\mathrm{exp}\{i(kx\omega t)\},p=p_0\mathrm{exp}\{i(kx\omega t)\},$$

where $`k`$ is the spatial wave number, $`\omega `$ the temporal frequency; $`u_0`$ and $`p_0`$ are respectively the amplitude of displacement and pressure. Upon substitution, it follows that:

$`p_0=(m\omega ^2ir\omega +s)u_0,`$
$`(k^2m+N)\omega ^2+i(k^2r+ϵ)\omega sk^2=0.`$ (2.7)
If the damping coefficients $`r=ϵ=0`$, then:

$$\omega =\pm \frac{ks_0^{1/2}}{(k^2m+N)^{1/2}},$$
(2.8)
showing that the system is dispersive, see , i.e, $`\omega =\omega (k)`$ and $`\omega ^{\prime \prime }(k)0`$. Waves of different wave length ($`2\pi /k`$) travel at different velocities. The dispersion relation (2.8) gives the group velocity:

$$\omega ^{}(k)=\pm s_0^{1/2}N(mk^2+N)^{3/2},$$
(2.9)
which decays to zero as $`k\mathrm{}`$. In other words, short waves ($`k`$ large) do not disperse as fast as long waves ($`k`$ small) .

When damping coefficients $`r`$, $`ϵ>0`$ ($`\beta =0`$), we have a damped dispersive system, (2.8) is modified into:

$$\omega =\frac{i(k^2r+ϵ)\pm \sqrt{4s_0k^2(k^2m+N)(k^2r+ϵ)^2}}{2(k^2m+N)},$$
(2.10)
with $`Im\{\omega \}=(k^2r+ϵ_s)/2(k^2m+N)<0`$, if the discriminant is nonnegative; otherwise:

$$Im\{\omega \}=\frac{(k^2r+ϵ_s)\pm \sqrt{4sk^2(k^2m+N)+(k^2r+ϵ_s)^2}}{2(k^2m+N)}<0.$$

In either case, we see decay of planar waves. The decay rate $`Im\{\omega \}`$ shows that if $`ϵ_s=0`$, the decay of long waves ($`|k|1`$) is very slow.  Hence the decay of long waves relies on $`ϵ_s`$.

We remark that the cochlear model in  also contains a so called stiffness coupling term in equation (2.2), so the elastic response becomes: $`p=mu_{tt}+r(x,t,u)u_t+s(x)uK(x)u_{xx},`$ for some positive function $`K(x)`$, rendering the highest spatial derivative order four in the model. Nevertheless, it is easy to check that the additional term does not change the dispersive nature of model solutions, neither the long wave dispersion and slow decay properties.

As we will see, the long wave dispersion and slow decay properties persist when coefficients become variables, and could cause problems for numerics. However, the generic dispersive effect of the time dependent solutions is absent in the special case when the temporal dependence of solution can be factored out. This is the case when nonlinearity is absent, i.e. $`r=r(x)`$. Let the input $`P_T(t)`$ be a single tone, or in complex variables $`P_T(t)=A\mathrm{exp}\{i\omega t\}+c.c`$. Writing solution as $`u=U(x)\mathrm{exp}\{i\omega t\}+c.c.`$, $`p=P(x)\mathrm{exp}\{i\omega t\}+c.c`$, we derive the following boundary value problem:

$`P(x)=(m\omega ^2+ir(x)\omega +s(x))U,`$ (2.11)
$`P_{xx}+(N\omega ^2ϵ(x)i\omega )U=0,`$ (2.12)
subject to the boundary conditions:

$$P_x(0)=T_MA,P(L)=0.$$
(2.13)
This is the way frequency domain calculation is done, turning a time dependent problem into a steady state problem on wave amplitudes, avoiding transients and the dispersion effect. However, when nonlinearity is present, it is in general impossible to factor out the time dependence exactly, especially when the input signal consists of multiple tones. To do an accurate computation directly on system (2.1)-(2.6), one has to proceed in the time domain, or solving the initial boundary value problem in time until transients die out and solutions approach a limiting time dependent state (quasiperiodic in general for multiple tones).

## 3 Numerical Method

In this section, we discuss our numerical method for solving system (2.1)-(2.6) and related numerical issues. Solutions of cochlear model with tonal inputs behave like traveling waves during early time, and eventually settle down to time periodic states of different periods at various BM locations according to the resonance information (frequency to place mapping) encoded in the function $`s(x)`$. A time domain computation will provide the stabilized limiting multi-frequency oscillating states after a long enough time integration. An ideal numerical method would be one that can reach the limiting states with accuracy, low cost and speed.

Let us first briefly review some existing methods. For a cochlear model with locally nonlinear damping function $`r`$ (van der Pol type), a method of line discretization is implemented in . Discretization in space is a Galerkin finite element method (of second order, with piecewise linear basis), and time stepping is explicit variable step 4th order Runge-Kutta, found to be superior to earlier explicit time stepping schemes such as Heun’s and modified Sielecki methods. A fully discrete explicit finite difference discretization, first order in time and second order in space, is carried out in . The discretization can be viewed as central differencing in space and forward Euler in time. The explicit nature of these methods put a severe stability constraint on the time step, requiring long computation time to reach stable quasi-periodic states that we are most interested in. Also it is known that one step explicit methods such as Runge-Kutta, Heun and Euler, are prone to accumulation of truncation errors for approximating dynamic objects like limit cycles, see , chapter 2.

Though implicit methods typically relax the stability constraints on time step and are better for computing steady states, they are also likely to be slow due to Newton iteration at each time step, and the non-banded complicated Jacobian matrix as a result of the nonlocal damping term. When $`r`$ is only Lipschitz continuous in $`u`$, lack of smoothness is another concern for the condition of the associated Jacobian matrix. This motivates us to look into semi-implicit methods that allow larger time steps and better stability properties than explicit methods, while being cheaper and faster than fully implicit schemes.

Let us introduce our 2nd order accurate semi-implicit discretization. We denote by $`u_j^n`$ ($`p_j^n`$) the numerical approximation of solution $`u(jh,nk)`$ ($`p(jh,nk)`$), where $`h`$ the grid size for $`x[0,L]`$, $`k`$ the time step. Let $`\delta _{\pm ,h}`$ be the forward/backward finite differencing operator in $`x`$, $`\delta _{\pm ,k}`$ the forward/backward finite differencing operator in $`t`$. Central first differencing operator in space is $`\delta _{0,h}=\frac{\delta _{+,h}+\delta _{,h}}{2}`$, and central second differencing operator in space is $`\delta _h^2=(\delta _{+,h}\delta _{,h})/h`$; with similar notations $`\delta _{0,k}`$, $`\delta _k^2`$ for central differencing operator in time. The semi-implicit method for system (2.1-2.6) is:

$`N\delta _k^2u_j^n={\displaystyle \frac{1}{4}}\delta _x^2(p_j^{n+1}+2p_j^n+p_j^{n1})ϵ_j\delta _{0,k}u_j^n,`$ (3.1)
$`p_j^{n+1}=m{\displaystyle \frac{2u_j^{n+1}5u_j^n+4u_j^{n1}u_j^{n2}}{k^2}}+r_j^n(u_t)_j^{n+1}+s_ju_j^{n+1},`$ (3.2)
where $`1jJ`$, $`n2`$, $`ϵ_j=ϵ_s(x_j)`$, and:

$$r_j^n=[r_{0,j}+\gamma h(\frac{1}{2}|u_1^{n+1}|e^{|x_jx_1|/\lambda }+\frac{1}{2}|u_J^{n+1}|e^{|x_jx_J|/\lambda }+\underset{l=2}{\overset{J1}{}}|u_l^{n+1}|e^{|x_jx_l|/\lambda })],$$
(3.3)
where $`(u_t)_j^{n+1}`$ is approximated at second order by:

$`(u_t)_j^{n+1}`$ $`=`$ $`(u_t)_j^{n1}+(2k)(u_{tt})_j^{n1}+O(k^2)`$ (3.4)
$`=`$ $`\delta _{0,k}u_j^{n1}+{\displaystyle \frac{2k}{N}}((p_{xx})_j^nϵ_j(u_t)_j^{n1})+O(k^2)`$
$`=`$ $`\delta _{0,k}u_j^{n1}+{\displaystyle \frac{2k}{N}}(\delta _h^2p_j^{n1}ϵ_j\delta _{0,k}u_j^{n1})+O(k^2),`$
and:

$`|u_l^{n+1}|`$ $`=`$ $`|u_l^n+k(u_t)_l^n+O(k^2)|`$ (3.5)
$`=`$ $`|u_l^n+(4u_l^{n1}u_l^{n2}3u_l^n)/(2)|+O(k^2)`$
$`=`$ $`|(4u_l^{n1}u_l^{n2}5u_l^n)/(2)|+O(k^2).`$
In (3.4), we have used the equation (2.1) once to lower time derivative of $`u`$ by one order. In (3.2), (3.3) and (3.5), we have also used one sided second order differencing to approximate $`u_{tt}`$, $`u_t`$, and trapezoidal rule to approximate the spatial integral. The variables at $`t=(n+1)k`$ are all linear. If follows from (3.2) that:

$$u_j^{n+1}=\frac{k^2p_j^{n+1}m(5u_j^n+4u_j^{n1}u_j^{n2})+kr_j^n(4u_j^nu_j^{n1})/2}{2m+3kr_j^n/2+k^2s_j}.$$
(3.6)
Substituting (3.6) into (3.1) to eliminate $`u^{n+1}`$, we get the following linear system of equations on $`p^{n+1}`$ ($`\lambda =k/h`$):

$`[(N+{\displaystyle \frac{ϵ_jk}{2}})k^2/(2m+{\displaystyle \frac{3kr_j^n}{2}}+k^2s_j){\displaystyle \frac{\lambda ^2}{2}}]p_j^{n+1}+{\displaystyle \frac{\lambda ^2p_{j+1}^{n+1}}{4}}+{\displaystyle \frac{\lambda ^2p_{j1}^{n+1}}{4}}`$ (3.7)
$`=`$ $`N(2u_j^n+u_j^{n1})+(N+{\displaystyle \frac{ϵ_jk}{2}}){\displaystyle \frac{m(5u_j^n+4u_j^{n1}u_j^{n2})+\frac{kr_j^n}{2}(4u_j^nu_j^{n1})}{(2m+\frac{3kr_j^n}{2}+k^2s_j)}}`$
$`{\displaystyle \frac{\lambda ^2}{4}}[2(p_{j+1}^n2p_j^n+p_{j1}^n)+p_{j+1}^{n1}2p_j^{n1}+p_{j1}^{n1}]{\displaystyle \frac{ϵ_jk}{2}}u_j^{n1},`$
which is solved by inverting a diagonally dominant tridiagonal matrix. Once $`p^{n+1}`$ is computed, $`u^{n+1}`$ is updated from (3.6).

The first two time steps are initiated as follows. Initial condition, (2.1), and (2.2) give:

$$p_{xx}\frac{N}{m}p=0,p_x(0)=f(0),p(L,0)=0,$$

whose solution is:

$$p(x,0)=\frac{f(0)}{\sqrt{N/m}(1+e^{2L\sqrt{N/m}})}(\mathrm{exp}\{(x2L)\sqrt{N/m}\}\mathrm{exp}\{\sqrt{N/m}x\}).$$
(3.8)
At $`t=k`$, we found

$$u(x,k)=k^2u_{tt}/2+O(k^3)=\frac{k^2}{2m}p(x,0)+O(k^3);$$

and similarly:

$$u_t(x,k)=kp(x,0)/m+k^2u_{ttt}(x,0)/2+O(k^3).$$

To find $`p(x,k)`$, denote $`q(x)=p_t(x,0)=mu_{ttt}(x,0)+r(x,0)u_{tt}(x,0)`$. It follows from (2.1) that $`q_{xx}Nu_{ttt}(x,0)=ϵ(x)u_{tt}(x,0)`$, implying

$$q_{xx}\frac{N}{m}q=ϵ(x)u_{tt}(x,0)\frac{N}{m}r(x,0)u_{tt}(x,0),$$

or a uniquely solvable two point boundary value problem on $`q`$:

$$q_{xx}\frac{N}{m}q=(Nr(x,0)/m^2ϵ(x))p(x,0),$$
(3.9)
with boundary condition: $`q_x(0)=f^{}(0)`$, $`q(L)=0`$.

It follows then from the definition of $`q`$ that:

$`u_{ttt}(x,0)`$ $`=`$ $`q(x)/mr(x,0)p(x,0)/m^2,`$
$`u_t(x,k)`$ $`=`$ $`kp(x,0)/m+{\displaystyle \frac{k^2}{2}}(q(x)/mr(x,0)p(x,0)/m^2).`$ (3.10)
Equations (2.1)-(2.2) at $`t=k`$ implies the uniquely solvable two point boundary value problem on $`p(x,k)`$:

$$p_{xx}\frac{N}{m}p=(ϵ(x)Nr(x,u))u_tNs(x)u/m,$$
(3.11)
with boundary condition: $`p_x(0)=f(k)`$, $`p(L)=0`$. Both (3.9) and (3.11) are solved with a standard second order discretization. We then will have computed $`(u,p)(x,k)`$. Similarly, with $`k`$ replaced by $`2k`$, we compute $`(u,p)(x,2k)`$.

The left boundary condition $`p_x(x,0)=f(t)`$ is discretized as: $`p_2^{n+1}p_0^{n+1}=2hf((n+1)k)`$, with second order accuracy. The right boundary condition $`p_{J+1}^{n+1}=0`$ is exact. The semi-implicit method is now completely defined, and it reduces to an unconditionally stable implicit method for standard second order wave equation when $`m=0`$, $`s`$ is constant, and damping is absent, , chapter 8.

## 4 Numerical Results and Multitone Interaction

In this section, we present numerical results based on our model and numerical methods discussed in the last two sections. First we list all parameters of our cochlear model and those of the associated middle ear model based on Guinan and Peake . We then show computed single tones and demonstrate the dispersion effects of long waves and the role of selective longitudinal fluid damping in our model. As a test of model tuning property, we also show the computed isodisplacement curves (at characteristic frequencies 500 Hz, 2 kHz, 4 kHz, 6 kHz, 10 kHz). Subsequently, we give numerical examples of both the high side and low side suppression of two tone interaction in model solutions consistent with Figures 10.2, 10.4, 10.7 in . The asymmetry of high and low side suppression is captured by the model. We then show numerical simulation of three tone interaction, presenting results on how the amplitudes of the two of the three tones change as we vary the amplitude and frequency of the third tone, qualitatively consistent with experimental data of Deng and Geisler  on responses of auditory neural fibers to input of three tones.

### 4.1 Model Parameters

The parameters for (2.2) are: cochlear length $`3.5cm`$; mass density $`0.05g/cm`$; $`\gamma =0.2`$; $`\lambda =0.2cm`$; $`N=16.67(dyne/cm^3)`$; $`r_0=0.001`$ $`g/(cm^2ms)`$.

The middle ear serves as a low pass filter from sound pressure at eardrum to the displacement of stapes, and as a high pass filter from the sound pressure at eardrum to the acceleration. We fit the data in  with the following gain factor $`a_M(\omega )`$ for each input $`e^{i\omega t}`$:

$$a_M(\omega )=30(\frac{1}{30}+0.0605\omega ^2((1\frac{\omega ^2}{\omega _m^2})^2+(2\xi _m\omega /\omega _m)^2)^{1/2}),$$
(4.1)
where $`\omega _m=4kHz`$, the middle ear characteristic frequency, and $`\xi _m=0.7`$ the middle ear damping ratio.

### 4.2 Dispersion Effect

We illustrate the dispersion property of time dependent solutions numerically in case of one tone input. In top plot of Fig. 1, we show the profiles of computed traveling wave at $`t=16`$ ms with and without the selective damping term, $`ϵ(x)u_t`$. For the run, $`h=0.01`$, $`k=0.01`$, and $`x_s=3`$, $`\beta =4`$. We see that for a tone input of $`4kHz`$, $`50dB`$, without selective damping, the BM displacement at $`t=16`$ ms contains a long wave which manages to pass through the characteristic location in the interval $`(1.5,2)`$ and persists near the exit $`x=L`$. This is due to the dispersive long waves with slow decay as we discussed in the last section. In contrast, when the selective damping term is present, the long wave is damped as it moves close to the right end point $`x=L`$. As the selective damping term is only effective near $`x=L`$, the solutions in the interior of the domain are essentially the same. In the bottom plot of Fig. 1, we show the BM displacement at $`t=32`$ ms, the dispersive long wave persists. In other runs (not shown here), we also observed growth in amplitude of long waves near $`x=L`$ with time at a linear rate. The appearance of long waves seems to be independent of numerical methods, rather their existence is intrinsic to the model which is dispersive.

### 4.3 Numerical Parameters

For all our runs reported below, $`h=k=0.01`$, $`ϵ=500`$, $`\beta =4`$, $`x_s=3`$. The $`x_s`$ and $`\beta `$ values should be properly increased ($`x_s`$ closer to $`x=3.5`$) if there are low frequencies below $`300Hz`$ in the input frequencies. Smaller $`h`$ and $`k`$ values have been used to test that numerical solutions do not vary much with further refinement. The time step $`k`$ should be made suitably smaller if input frequencies are as large as $`16kHz`$ for more resolution.

### 4.4 Isodisplacement Curves

For each input pressure wave of frequency $`\omega `$, there is a unique location $`x_\omega [0,L]`$ where the maximum of BM displacement (peak in absolute value) is located. The $`x_\omega `$ is called characteristic place. Conversely, for each location $`x[0,L]`$ there is a frequency, called characteristic frequency (CF) and denoted by $`CF(x)`$, such that the BM response amplitude attains the maximum at $`x`$. In Fig. 1, we see that $`x_\omega (1.5,2)`$ for $`\omega =4`$ kHz. The higher the amplitude of the input wave (denoted by $`A_{in}`$), the higher the BM response peak. If the level of BM response is set at a fixed level, such as 1 nm, and a BM location $`x`$ is given, we look for an input value $`A_{in}`$ so that $`|u|(x)=1`$ nm. Here $`|u|`$ is the steady state BM response for input wave $`A_{in}e^{i\omega t}`$. The plot of $`A_{in}`$ as a function of $`x`$ gives the so called isodisplacement curve. The profile $`|u|(x)`$ is asymmetric about the peak, and decays rapidly to zero beyond $`x_\omega `$. So it takes higher $`A_{in}`$ to stimulate a point $`x>x_\omega `$ to 1 nm than to the left. It takes least $`A_{in}`$ to stimulate $`|u|(x_\omega )`$ to 1 nm. So the isodisplacement curve has a minimum at $`x_\omega `$, rises sharply at $`x>x_\omega `$ and gradually at $`x<x_\omega `$.

In Fig. 2, we show model isodisplacement curves for five characteristic frequencies 500 Hz, 2 kHz, 4 kHz, 6 kHz, 10 kHz. The threshold displacement is 1 nm. The curves showed the trend of sharper tips in high frequencies, wider tips for lower frequencies, and the asymmetry about the tips. The plot is comparable to the one in . Isodisplacement curves are related to frequency selectivity of hearing.

### 4.5 Two and Three Tone Interactions

Two tone interaction is well-documented in the literature, here we illustrate that our model gives qualitatively the same results as shown in related figures in . We shall use BM displacement to detect the tonal interaction. Fig. 3 shows that the amplitude of a 5 kHz 70 dB tone drops (by nearly 50%) after the second tone of 2.4 kHz and 70dB is introduced, this reveals the so called low side suppression.

In Fig. 4, we see that the 5 kHz tone at 50 dB (solid line) is suppressed (by about 30 %) after interacting with 6.7 kHz tone at 80 dB (plus-line), the so called high side suppression. The high side suppressor tone at 6.7 kHz has to be much higher in amplitude than a low side suppressor tone whose frequency is at least as far away from that of the suppressed tone. This shows the asymmetry of low and high side suppression as in  among others.

In , Fig. 7a showed cat auditory neural rate response at a fixed characteristic frequency (CF) demonstrating suppression effect for two tone $`(F_1,F_2)`$ input, where $`F_1`$ at fixed intensity takes three frequency values from below CF to above CF, for various values of the frequency $`F_2`$ at 70 dB. In Fig. 5, we showed a qualitatively similar plot computed with our model, where CF=1 kHz, 50 dB $`F_1`$ takes three frequency values (0.9, 1, 1.2) kHz, 70 dB $`F_2`$ increases its frequencies as 1.3 kHz, 1.4 kHz, 1.5 kHz, 1.6 kHz, 1.7 kHz, 1.8 kHz.

The three tone interaction is much less documented in the literature, here we shall compare our model results qualitatively with the experimental data on auditory neural response to multi-tones in . Two tones are fixed at 4 kHz, 4.4 kHz, both at 50 dB. The third tone varies in frequency and in amplitude. In Fig. 6, we plot the BM response of two tones (4 kHz, plus line; 4.4 kHz, dotted line) as the third tone (line) increases its amplitude from 33 dB to 73 dB for three values of frequencies, 3.6 kHz, 3.8 kHz, 4.2 kHz. We see that the two tone response decreases for increasing amplitude of third tone; and that the suppression effect on the 4kHz tone is larger at 3.8 kHz (low side) than at 4.2 kHz (high side) for the same amplitude of third tone. The effect of third tone at 3.6 kHz is even larger than that at 4.2 kHz. In Fig. 7, we plot similarly the response of two tones as the third tone takes on frequencies 4.3 kHz, 4.6 kHz, 4.8 kHz. The masking effect of third tone remains, except that the second tone at 4.4 kHz shows more difference as we observe the crossing of two tone curves when the frequency of third tone goes through that of the second tone (4.4 kHz). The third tone is more efficient in masking the second tone from the low side (4.3 kHz) than from the high side (4.6 kHz). The suppressing effect on the second tone is even weaker at 4.8 kHz. These features agree with those of Fig. 3 of Deng and Geisler , e.g. compare frame 4 and frame 6 there. Strictly speaking the neural data in  are synchrony masking responses, and BM responses we computed correspond to rate masking responses. However, their data are relevant for qualitative comparison, as there is close (positive) correlation between the rate and phase behavior, see .

## 5 Conclusions

We discussed the dispersive property of the cochlear models for time domain computation, and slow decay of long waves. A selective longitudinal damping term has been introduced in the nonlinear nonlocal model studied previously . The new model is computed using a semi-implicit second order finite difference method in the time domain. We presented numerical results on isodisplacement curves, two tone suppression, and three tone interaction in qualitative agreement with earlier findings of Geisler  and experimental auditory neural data Deng and Geisler . These encouraging results prompt one to perform further study of masking mechanism in more complex tones and speech input based on our cochlear model.

## 6 Acknowledgments

J.X would like to thank Prof. C. D. Geisler for helpful email communications on two tone suppression and cochlear modeling. He wishes to thank Prof. J. Keller for his critical reading of the first manuscript and his constructive comments, and Prof. S. Osher for his encouragement. The authors wish to thank the Institute for Mathematics and its Applications (IMA) at the University of Minnesota for hosting the workshop on speech processing in Sept of 2000 and for providing an interactive scientific environment, where the initial collaboration began. The authors also thank Mr. M. D. LaMar for his interest and careful reading of the paper.