format: text

text:
# Localization of Denaturation Bubbles in Random DNA Sequences

## I Introduction

Localized opening of double-stranded DNA is essential in a number of cellular processes such as the initiations of gene transcription and DNA replication Alberts . While thermal denaturation is highly unlikely under physiological conditions, in vitro experiments show that local denaturation can be readily induced by under-winding the DNA double helix by an amount that is physiologically reasonable Kowalski , Cozzarelli , Bensimon . The basic physical effect is simple: An under-wound double helix suffers a reduction in the binding free energy Marko , Cocco , Strick00 . Local openings of the double helix (referred to as “denaturation bubbles”) relieve the twist experienced by the remainder of the double helix and is thus energetically favorable. The denaturation bubbles may be recruited to a specific location of the genome by a designed (e.g., AT-rich) sequence, since AT pairs bind more weakly than GC pairs Benham . On the other hand, entropic effect which favors bubble delocalization is non-negligible for long sequences. Also significant is the kinetic trapping of the bubbles due to statistical agglomeration of AT-rich segments in long heterogenous sequences.

To gain some quantitative understanding on the competing effects of entropy and sequence heterogeneity, we characterize in this study the thermodynamic and dynamic properties of denaturation bubbles in a long, stretched random DNA sequence with no special sequence design. Previously, there have been a number of experimental and theoretical studies Cule , tang01 , Heslot , Lubensky  on the effect of sequence heterogeneity on DNA melting and unzipping transitions. Our study is along this general direction. The specific behaviors exhibited by the denaturation bubbles are rather complex and are typical of those observed in systems dominated by quenched disorders glass : The bubbles are localized upon increase of the applied torque beyond a certain threshold. In the localized regime, their dynamics exhibits “aging” aging1 , aging  and is sub-diffusive with continuously varying exponents.

Interestingly, twist-induced denaturation presents a rare physical example of the celebrated random-energy model of disordered system rem . Consequently, detailed analysis of both the thermodynamis and dynamic properties can be made by applying the well-developed theory of disordered systems glass , together with exact results from large-deviation theory familiar in the related sequence alignment problem KarlinAltschul , YuHwa . We will draw upon detailed experimental knowledge of thermal denaturation SantaLucia , SantaLucia1 , Blake  throughout the analysis, and make our results quantitative whenever possible.

## II Thermodynamics

Let us consider the application of a torque which under-winds a long, stretched111A modest stretching force is needed to prevent the applied torque from being absorbed by super-coiling; see e.g., Ref. Cocco . piece of double-stranded DNA. We are interested in the regime where the applied torque $`𝒯`$ is below the threshold $`𝒯_\mathrm{d}`$ for bulk denaturation, but sufficiently strong so that denaturation bubbles do appear in the system. Due to the highly cooperative nature of the denaturation process, the typical distance $`N_\times `$ between the large bubbles is large, in which case treating the bubbles as a dilute gas of particles is appropriate. Our strategy will be to first characterize analytically the thermodynamic behavior of a single bubble, and then use this knowledge to determine the length scale $`N_\times `$ and the many-bubble states for $`NN_\times `$. We will find that $`N_\times >O(10^3)\mathrm{𝚋𝚙}`$ as long as we are not very close to the threshold $`𝒯_\mathrm{d}`$, so that the dilute gas approximation is reasonable for a large range of parameters.

### II.1 The Single-bubble Model

Consider a denaturation bubble confined in a DNA double helix between two complementary DNA strands of $`N`$ bases each. The double-strand is denoted by the base sequence $`b_1b_2\mathrm{}b_N`$ \[with $`b_k\{𝙰,𝙲,𝙶,𝚃\}`$\] of one of the strands, ordered from the 5’- to 3’- end. To simplify the notation, we assume that the two ends of the helix is sealed, so that the bubble is always contained in the segment $`b_1\mathrm{}b_N`$. Let the index of the first and last open pairs of the open bubble be $`m`$ and $`n`$, with $`1mnN`$. We denote the total free energy of the bubble (defined with respect to the helical state) by $`\mathrm{\Delta }G_L(m)`$, where $`Lnm+1`$ is the number of open bases and referred to as the bubble length. Then the partition function of the single-bubble system is given by

$$Z(N)=\underset{L=1}{\overset{N}{}}\underset{m=1}{\overset{NL+1}{}}e^{\beta \mathrm{\Delta }G_L(m)},$$
(1)
where $`\beta ^1k_BT0.62\mathrm{𝚔𝚌𝚊𝚕}/\mathrm{𝚖𝚘𝚕𝚎}`$ at $`37^{}𝙲`$.

In the absence of the external torque, the bubble energy $`\mathrm{\Delta }G_L(m)`$ has two components. First, there is the loss of stacking energy $`\delta G_{b,b^{}}`$ between two successive bases $`b`$ and $`b^{}`$. These stacking energies are in the range $`0.52.5k_BT`$’s at $`37^{}𝙲`$, with the AT stacks weaker than the GC stacks. Their values have been measured carefully SantaLucia , SantaLucia1 , Blake . Second, assuming that there is no secondary pairing between bases in the bubble so that the open configuration can be regarded as a polymer loop, there is a well-known polymeric loop entropy cost

$$\gamma _L=\gamma _1+\alpha k_BT\mathrm{ln}L$$
(2)
for a bubble of length $`L`$, with $`\alpha 1.8`$ MEFisher  for a linearly extended 222The value of $`\alpha `$ may well be different for unstretched DNA chain and hence relevant for the thermal denaturation of DNA Mukamel , Orland . However, as we show below, essential features of denaturation do not hinge on the precise value of $`\alpha `$. DNA chain. The bubble initiation cost $`\gamma _1`$ depends on the base composition at opening and closing ends, ionic strength etc., and generally lies333The initiation cost for DNA bubbles are extracted from the webserver: http://www.bioinfo.rpi.edu/applications/mfold/ (M. Zuker, private communication). See also Ref. Rouzina  for an alternative source. in the range $`35k_BT`$. For relevant bubble sizes of few tens of bases in length (see below), the total entropic cost is $`\gamma _L=812k_BT`$. This large cost justifies the single bubble approximation (at least to the length scale $`N_\times e^{\beta \gamma _L}5\times 10^3\mathrm{𝚋𝚙}`$), and contributes significantly to the sharpness of the observed thermal denaturation transition Blake .

An applied negative torque $`𝒯`$ reduces the thermodynamic stability of the helical state relative to the denatured one by an amount equal to the work done to unwind the helix. This effect is simply modeled here by a linear decrease in the stacking energy in the relevant parameter range Cocco , i.e., $`\delta G_{b,b^{}}\delta G_{b,b^{}}\theta _0𝒯`$, where $`\theta _0=2\pi /10.35`$ is the twist angle per base of the double helix. Putting the above together, we have

$$\mathrm{\Delta }G_L(m)=\gamma _L+\underset{k=m}{\overset{m+L}{}}\delta G_{b_{k1},b_k}\theta _0𝒯L$$
(3)
as the single bubble energy, which can be computed once the DNA sequence $`b_1\mathrm{}b_N`$ is given. Note that while Eq. (3) is formulated specifically for twist-induced denaturation, the general form can be used to describe a number of destabilizing effects, e.g., due to changes in temperature, ionic concentration, etc.

### II.2 Sequence Heterogeneity

As the torque $`𝒯`$ increases from zero towards the denaturation point $`𝒯_\mathrm{d}`$, denaturation bubbles appear in the double strand and grow in size. We wish to know whether the bubbles are free to diffuse along the double strand, or are they localized in the high AT regions of the DNA where binding is the weakest. For simplicity, we will characterize the typical behavior of an ensemble of random (i.e., independent and identically distributed) sequences described by the single-nucleotide frequencies $`p_b`$, although our approach and qualitative findings are also applicable to sequences with short-range correlations.

For a given sequence of bases, the partition function $`Z`$ can of course be efficiently evaluated numerically (including all the multiple-bubble states) by using available programs such as MELTSIM Blake . All thermodynamic quantities can subsequently be evaluated from the free energy $`F=k_BT\mathrm{ln}Z`$. To obtain the typical behavior of the ensemble, we ideally want to compute the ensemble average of the free energy, $`\overline{F}k_BT\overline{\mathrm{ln}Z}`$. \[We use the over-line to denote average over the ensemble of random sequences, i.e., $`\overline{X}_{b_1,\mathrm{},b_N}X_{b_1,\mathrm{},b_N}_{k=1}^Np_{b_k}`$; this is also known as the “disorder average”.\] Computing $`\overline{F}`$ numerically however will require explicit generation of a large number of random sequences and can be very time consuming for large $`N`$’s. Fortunately we can apply a large body of knowledge accumulated from the statistical mechanics of random systems glass  and provide a detailed characterization of the typical behavior of our system without the need of exhaustive simulation. To introduce notation and concepts in this approach, we examine first the simplified problem of a single bubble with a fixed length.

### II.3 Bubble with Fixed Length

Let us consider a bubble with a fixed length $`L`$ (with $`1LN`$) embedded in a long, random sequence $`b_1\mathrm{}b_N`$. The partition function reads

$$𝒵_L=\underset{m=1}{\overset{NL+1}{}}\mathrm{exp}[\beta \mathrm{\Delta }G_L(m)],$$
(4)
where the scripted variables refer to properties of the fixed length bubble. For a random sequence, the energies of the different states labeled by $`m`$ are uncorrelated with each other beyond the distance $`L`$. Such systems belong 444The correlation in $`\mathrm{\Delta }G`$ between neighboring states is only a minor complication because it is short-ranged and can be transformed away by coarse graining. to the class of “Random Energy Model” (REM) and was solved exactly in the 80’s by Derrida rem  for a Gaussian distribution of $`\mathrm{\Delta }G`$’s. Discrete distribution of $`\mathrm{\Delta }G`$’s was studied in the closely-related system involving protein-DNA interaction Gerland . Below, we will briefly review the salient properties of REM using the present example.

The REM has a “high-temperature” phase where many (of order $`N`$) bubble configurations contribute significantly to the partition sum, and a “low-temperature” phase dominated by only one or a few lowest energy states. It follows that, in the former case, the bubble is delocalized and can freely diffuse along the sequence, while in the latter case, the bubble is localized to the lowest energy position. Transition between the two phases is driven by competition between the energetic (variation in $`\mathrm{\Delta }G`$) and entropic ($`\mathrm{ln}N`$) effects. In the present problem, the magnitude of terms in the partition sum (4) can be tuned not only by varying the temperature, but also by varying the bubble size $`L`$. Hence, at a fixed $`\beta `$, whether a bubble is free or localized depends both on the bubble size $`L`$ and the sequence length $`N`$.

An interesting property of the REM is that, in the “high-temperature” phase, $`𝒵_L/N`$ tends to a finite limit given by the annealed average $`\overline{𝒵}_L/N`$ as $`N\mathrm{}`$, independent of the particular realization of the random sequence. This allows us to replace the average free energy $`\overline{}k_BT\overline{\mathrm{ln}𝒵}`$ by its annealed approximation $`\stackrel{~}{}k_BT\mathrm{ln}\overline{𝒵}`$, which is much easier to calculate. \[We will use the tilde to denote all quantities computed in the annealed approximation.\] Introducing a $`4\times 4`$ matrix $`𝗠(\beta )`$ with components $`𝖬_{b,b^{}}=\sqrt{p_bp_b^{}}\mathrm{exp}[\beta \delta G_{b,b^{}}]`$, and let the largest eigenvalue of $`𝗠(\beta )`$ be $`\mathrm{\Lambda }(\beta )`$, then the disorder average of terms in (4) can be written as

$$\overline{e^{\beta _{k=m}^{m+L1}\delta G_{b_{k1},b_k}}}=\mathrm{Tr}𝗠^L(\beta )\stackrel{L1}{=}\mathrm{\Lambda }^L(\beta ).$$
(5)
It is convenient to introduce the quantity

$$f(\beta )\beta ^1\mathrm{ln}\mathrm{\Lambda }(\beta ),$$
(6)
with which we have (for $`NL`$)

$$\overline{𝒵}_L=\frac{Ne^{\beta \gamma _1}}{L^\alpha }\left[e^{\beta [f(\beta )\theta _0𝒯]}\right]^L.$$
(7)
Hence, in the delocalized phase, $`\overline{}\stackrel{~}{}=k_BT\mathrm{ln}N+L[f(\beta )\theta _0𝒯]+\gamma _L`$.

The annealed entropy can be calculated from $`\stackrel{~}{}`$, with555To focus on the positional entropy, we did not include here the contribution due to loop entropy, i.e., we treated $`\gamma _L`$ as an energy term despite its entropic origin.

$$\stackrel{~}{𝒮}=k_B^1\frac{\stackrel{~}{}}{T}=\mathrm{ln}N\beta [f(\beta )\epsilon (\beta )]L,$$
(8)
where $`\epsilon (\beta )\frac{}{\beta }\mathrm{ln}\mathrm{\Lambda }.`$ It will also be useful to introduce the relative entropy per base for the fixed-length bubble,

$$(\beta )[\mathrm{ln}N\stackrel{~}{𝒮}]/L=\beta [f(\beta )\epsilon (\beta )].$$
(9)
Note that being the difference between $`f`$ and $`\epsilon `$, the quantity $``$ is a measure of the intrinsic variation in the binding energies $`\delta G`$’s for a random sequence with nucleotide frequency $`p_b`$, and is independent of the average binding energy $`\overline{\delta G}`$ which external environment such as the temperature or solvent conditions most directly affect.

Derrida’s solution of REM shows that the annealed entropy $`\stackrel{~}{𝒮}`$ vanishes at the transition to the “low-temperature” phase, beyond which the annealed approximation is no longer applicable. Using Eqs. (8) and (9), we can write the condition for phase transition as

$$L_{\mathrm{loc}}=\mathrm{ln}N/(\beta ),$$
(10)
which gives the minimal bubble size for localization at a given $`N`$. With the values of $`\delta G`$’s obtained from Ref. SantaLucia  at $`[\mathrm{Na}^+]=1𝙼`$ and $`37^{}𝙲`$, and assuming an equal nucleotide distribution (i.e., $`p_b=1/4`$ for all bases), we find $`f1.83k_BT`$, $`\epsilon 1.50k_BT`$, so that $`0.33`$ and $`L_{\mathrm{loc}}20\mathrm{𝚋𝚙}`$ for $`N10^3\mathrm{𝚋𝚙}`$. From Eq. (10), it is clear that as $`N\mathrm{}`$ any fixed-length bubble remains delocalized.

### II.4 Bubble without Length Constraint

The full partition function $`Z`$ is obtained simply by summing $`𝒵_L`$ for different $`L`$’s. We will again approach the problem by first applying the annealed approximation and then determining where it breaks down.

1. Annealed approximation: The annealed partition function $`\overline{Z}(N)_{L=1}^N\overline{𝒵}_L`$ has a transition at $`𝒯_\mathrm{a}=f(\beta )/\theta _0`$, where the exponential factor in (7) reaches one: The sum over $`L`$ is finite and $`\overline{Z}N`$ only if $`𝒯<𝒯_\mathrm{a}`$. In this regime, the annealed free energy is simply $`\stackrel{~}{F}k_BT\mathrm{ln}N+\gamma _1.`$ The annealed energy $`\stackrel{~}{E}\frac{}{\beta }\mathrm{ln}\overline{Z}(\beta )`$ is also readily computed; it can be expressed as $`\stackrel{~}{E}=\gamma _1+[\epsilon (\beta )\theta _0𝒯]\stackrel{~}{L}`$ where $`\stackrel{~}{L}(𝒯)_{L=1}^N\mathrm{}L\overline{𝒵}_L/\overline{Z}`$ is the average bubble length in the annealed approximation. As $`𝒯`$ approaches $`𝒯_\mathrm{a}`$, $`\stackrel{~}{L}(𝒯)`$ diverges, and the annealed entropy

$$\stackrel{~}{S}=\mathrm{ln}N\beta [\theta _0𝒯\epsilon (\beta )]\stackrel{~}{L}$$
(11)
becomes negative.

In the limit $`N\mathrm{}`$, the anealed free energy $`\stackrel{~}{F}`$ is actually identical to $`\overline{F}`$ for all $`𝒯𝒯_\mathrm{a}`$. This can be seen from the inequalities $`\overline{\mathrm{ln}𝒵_{L=1}}\overline{\mathrm{ln}Z}\mathrm{ln}\overline{Z},`$ and $`𝒵_{L=1}>N\mathrm{min}\{\mathrm{exp}[\beta \mathrm{\Delta }G_1]\}`$. Since both the lower and upper bounds grow as $`\mathrm{ln}N`$,

$$\overline{F}k_BT\overline{\mathrm{ln}Z}k_BT\mathrm{ln}N$$
(12)
for all $`𝒯𝒯_\mathrm{a}`$.

2. Ground-state properties: To find the ground-state of the unconstrained bubble in a long random sequence, we need to study the statistics of stretches of exceptionally high AT-content. If we neglect the polymeric contribution $`\gamma _L`$ to the bubble energy \[to be justified shortly\], then the ground-state energy $`E^{}`$ expected in a sequence of length $`N`$ can be computed exactly from large-deviation theory KarlinAltschul , Dembo , with

$$\overline{E^{}}(N)\lambda ^1\mathrm{ln}N.$$
(13)
The constant $`\lambda `$ in Eq. (13) can be expressed as the unique positive root of the equation

$$f(\lambda )=\theta _0𝒯,$$
(14)
where $`f`$ is defined by the $`\delta G`$’s through Eqs. (5) and (6). Note that, at $`𝒯=𝒯_\mathrm{a}`$, Eq. (14) is satisfied with $`\lambda =\beta `$. In this case, (13) coincides with (12).

The length of the minimal energy bubble is also known from large-deviation theory Dembo , YuHwa , with

$$L^{}(𝒯)=\mathrm{ln}N/H^{}(𝒯),$$
(15)
where the relative entropy $`H^{}`$ is given exactly by

$$H^{}(𝒯)=\lambda (𝒯)[\theta _0𝒯\epsilon (\lambda )].$$
(16)
From the logarithmic dependence of the bubble length $`L^{}`$ on $`N`$, it is clear that the corresponding polymeric contribution $`\gamma _L^{}\mathrm{ln}(\mathrm{ln}N)`$ can indeed be treated as a constant shift of bubble energy.

3. Phase transitions: Based on the above discussion, a phase transition can be formally established in the limit $`N\mathrm{}`$. This is seen by comparing the expressions (12) and (13). For $`𝒯>𝒯_\mathrm{a}`$, solution to (14) satisfies $`\lambda <\beta `$. Consequently, (12) must break down there, yielding a phase transition at $`𝒯=𝒯_\mathrm{a}`$. Since $`\overline{F}\overline{E^{}}`$ in general (e.g., for all $`𝒯>𝒯_\mathrm{a}`$), and at the phase transition point $`𝒯=𝒯_\mathrm{a}`$ the equality $`\overline{F}=\overline{E^{}}`$ already holds, i.e., the ground-state already dominates, then we must have the ground-state dominating throughout the localized phase. This is exactly the behavior of the random energy model rem .

A physical understanding of the transition can be obtained by examining the importance of the ground-state contribution $`\mathrm{exp}(\beta E^{})N^{\beta /\lambda }`$ to the partition sum $`Z`$ as the applied twist $`𝒯`$ is varied. For $`𝒯<𝒯_\mathrm{a}`$, the ratio $`\beta /\lambda (𝒯)`$ is less than one. In this case, the energy gain $`\overline{E^{}}(N)`$ of placing the bubble at the site of the lowest energy is insufficient to overcome the entropy $`\mathrm{ln}N`$ of placing the bubble in different positions, hence the bubble is typically small and delocalized. When $`𝒯`$ exceeds $`𝒯_\mathrm{a}`$, the ground-state contribution grows faster than $`N`$, signalling dominance of one or a few low energy states where the bubble typically resides. The transition is thus identified as the localization transition of the bubble at $`𝒯_{\mathrm{loc}}=𝒯_\mathrm{a}`$.

The onset of the zero entropy point can be obtained from Eq. (11) and written as

$$\mathrm{ln}N=H(\beta ,𝒯)\stackrel{~}{L}(\beta ,𝒯),$$
(17)
where

$$H(\beta ,𝒯)=\beta [\theta _0𝒯\epsilon (\beta )]$$
(18)
is the relative entropy of the unconstrained bubble. These equations are analogous to the expressions (15) and (16) for the ground-state bubble. In fact, both $`L^{}(𝒯)`$ and $`H^{}(𝒯)`$ are reproduced through the substitution $`\beta \lambda (𝒯)`$, e.g., $`L^{}(𝒯)=\stackrel{~}{L}(\lambda (𝒯),𝒯)`$. This turns out to be true also for other thermodynamic variables. Thus the localized phase at different $`𝒯`$’s can be viewed as the phase transition points of systems with different effective temperatures $`\lambda ^1(𝒯)`$; this will be clearly manifested in the bubble dynamics discussed below.

Next we observe that since $`H^{}\lambda `$ \[see Eq. (16)\], the bubble length diverges (or approaches $`N`$) as $`\lambda 0`$. This defines the point of bulk denaturation666Note however that the helical segments separating adjacent bubbles can be stable even beyond $`𝒯_\mathrm{d}`$, so that complete separation of the two strands takes place at $`𝒯>𝒯_\mathrm{d}`$.  $`𝒯_\mathrm{d}`$, i.e.,

$$\theta _0𝒯_\mathrm{d}\underset{\lambda 0}{lim}f(\lambda )=\overline{\delta G},$$
(19)
where the second equality is obtained from manipulating Eqs. (5) and (6). Using $`\overline{\delta G}1.40k_BT`$ \[derived from the $`\delta G`$’s in Ref. SantaLucia \], we find $`𝒯_\mathrm{d}10\mathrm{pN}\mathrm{nm}`$. The dependence of $`\lambda `$ on $`𝒯`$ close to $`𝒯_\mathrm{d}`$ can be obtained from the expansion

$$f(\lambda )=\overline{\delta G}\frac{\lambda }{2}\mathrm{var}(\delta G)+O(\lambda ^2).$$
(20)
Inverting the above for $`\lambda `$ using (14) and (19), we find

$$\lambda (𝒯)=\frac{2\theta _0}{\beta \mathrm{var}(\delta G)}\left(𝒯_\mathrm{d}𝒯\right)+O\left(𝒯_\mathrm{d}𝒯\right)^2.$$
(21)
It turns out that the term linear in $`𝒯_\mathrm{d}𝒯`$ in (21) already gives a very good approximation (to within $`1\%`$) of $`\lambda `$ throughout the localized phase where $`\lambda /\beta <1`$. The localization transition point $`𝒯_{\mathrm{loc}}`$ can be thus obtained by solving Eq. (21) with $`\lambda (𝒯_{\mathrm{loc}})=\beta `$. Using $`\beta ^2\mathrm{var}(\delta G)0.565`$ \[derived from Ref. SantaLucia \], we find $`𝒯_\mathrm{d}𝒯_{\mathrm{loc}}2\mathrm{pN}\mathrm{nm}`$. Unlike the value of $`𝒯_\mathrm{d}`$ which is derived from the average stacking energy $`\overline{\delta G}`$ and hence sensitive to temperature, ionic strength, etc., the difference $`𝒯_\mathrm{d}𝒯_{\mathrm{loc}}`$ is set by the variance of $`\delta G_{b,b^{}}`$ and should be much less sensitive to experimental conditions. The same is expected for the relative entropy, which has the form

$$H^{}(𝒯)2\theta _0^2(𝒯_\mathrm{d}𝒯)^2/\mathrm{var}(\delta G)$$
(22)
throughout the localized phase.

### II.5 Multiple bubbles

The localization transition discussed above occurs only as $`N\mathrm{}`$. However for large $`N`$, the single-bubble approximation will break down regardless of the large (but finite) bubble cost $`\gamma _L`$. When multiple bubbles are localized, each bubble is effectively in a finite-length system, thereby blurring the localization transition.

We first analyze the delocalized phase for which the annealed approximation is valid. Once multiple bubbles are allowed in the system, we expect a broad range of bubble lengths, as described by the distribution (7). Qualitatively, we expect only the largest bubbles, of size $`\stackrel{~}{L}(𝒯)`$ to be localized as $`𝒯𝒯_{\mathrm{loc}}`$, while the smaller ones remain delocalized. We shall thus focus on these large bubbles: It is the average separation distance $`N_\times `$ between these large bubbles that sets the effective system size of the single-bubble localization problem.

The Boltzmann weight of one such large bubble in a sequence of length $`N\stackrel{~}{L}`$ is $`\stackrel{~}{W}(N)=e^{\beta \gamma _1}N/\stackrel{~}{L}^\alpha `$ in the vicinity of the localization transition. Setting $`\stackrel{~}{W}(N)=1`$ yields the typical spacing between the large bubbles on the delocalized side,

$$\stackrel{~}{N}_\times e^{\beta \gamma _1}\stackrel{~}{L}^\alpha (𝒯).$$
(23)
Note that for bubbles of size $`10\mathrm{𝚋𝚙}`$’s, the crossover length is already of the order of $`10^3\mathrm{𝚋𝚙}`$’s. A similar estimate can be made on the localized phase using the exact expression KA93  for the lowest energy for multiple bubbles. We find

$$N_\times ^{}\left[e^{\beta \gamma _1}(L^{})^\alpha \right]^{\lambda (𝒯)/\beta }.$$
(24)
as the average distance between large bubbles of size $`L^{}`$.

For $`NN_\times `$, the system consists of $`N/N_\times `$ effective number of single-bubble subsystems, each of length $`N_\times `$. At the localization “transition” of an infinite system then, we have $`\mathrm{ln}\stackrel{~}{N}_\times =H(\beta ,𝒯_{\mathrm{loc}})\stackrel{~}{L}(𝒯_{\mathrm{loc}})`$ \[see Eq. (17)\]. Together with Eq. (23) \[or (24) with $`\lambda =\beta `$\], we find $`\stackrel{~}{L}(𝒯_{\mathrm{loc}})25\mathrm{𝚋𝚙}`$ at the onset of localization (using $`\gamma _13k_BT`$ and $`H0.33`$), with the crossover length $`\stackrel{~}{N}_\times 6500\mathrm{𝚋𝚙}`$. Thus we expect there to be typically one bubble of $`25\mathrm{𝚋𝚙}`$ in a random DNA double-strand of length $`6500\mathrm{𝚋𝚙}`$’s at the localization transtion.

## III Bubble Dynamics

The localization of bubbles is reflected ultimately in their slow dynamics. We expect bubbles to diffuse freely along the DNA double-helix in the delocalized phase, but become trapped in low energy positions in the localized phase. Details of the bubble movement in the latter case, however, can be rather complicated with nontrivial memory (or “aging”) effects typical of glassy states aging1 , aging  as will be described below.

### III.1 Model

For simplicity, we will restrict ourselves to the description of the movement of a single-bubble over its lifetime, which can be rather long in the localized phase. For reasons discussed above, interaction with other bubbles can be neglected when the bubble displacement is within a distance of order $`N_\times 10^3\mathrm{𝚋𝚙}`$. We will also neglect the polymeric loop entropy $`\gamma _L`$ which provides essentially a constant shift to the bubble energy as shown in the single bubble section.

In addition to the drift and breathing motion, a bubble may also shrink to zero size and disappear from the system. To our knowledge, the time scale involved for the spontaneous collapse of a bubble, particularly under an applied twist, has not been documented so far. Zipping the bubble requires not only pairing of the bases in the open segment, but also rewinding of the helix against the applied undertwist, both of which contribute to the energy barrier to the no-bubble state. This suggests a long lifetime for a bubble, which can be enforced by setting a lower bound (e.g., $`10\mathrm{𝚋𝚙}`$) in the allowed bubble length. However, as we will see, the long-time behavior of bubble dynamics is determined crucially by the occurrence of the large bubble states, and insensitive to the value of the lower bound on $`L`$, as long as the $`L=0`$ state is excluded. Once accurate estimates of bubble lifetime become available, one may supplement the discussion below with such a cutoff.

### III.2 Scaling Theory

Equation (13) gives the lowest energy of an unconstrained bubble in a sequence of length $`N`$, while a bubble with its position (but not size) fixed typically has an energy of the order $`\lambda ^1`$ for $`𝒯<𝒯_{\mathrm{loc}}`$. For small $`\lambda `$, the energy variation $`\mathrm{\Delta }E(N)\lambda ^1\mathrm{ln}N`$ is large, hence the bubble dynamics is dominated by the thermal escape from the deepest trap. The escape time is thus $`t_\mathrm{e}(N)e^{\beta \mathrm{\Delta }E(N)}N^{\beta /\lambda }`$, i.e., the dynamics is sub-diffusive deep in the localized phase (where $`\beta \lambda `$).

To investigate the dynamical behavior in more detail, especially close to the localization transition where $`\lambda \beta `$, we need to include also the random motion of the bubble along the double strand. Towards this end, it is useful to describe the bubble dynamics as a single point moving in the two-dimensional space spanned by the bubble’s only two degrees of freedom, its instantaneous length $`L`$ and the position of one of its ends, say $`m`$. The statistics of the two-dimensional energy landscape $`\mathrm{\Delta }G_L(m)`$ is well-characterized by the large-deviation theory YuHwa . It consists of a number of valleys, whose depths (denoted by $`\mathrm{\Delta }\widehat{G}`$’s) are given by the Poisson distribution

$$𝒟\left(\left|\mathrm{\Delta }\widehat{G}\right|\right)e^{\lambda \left|\mathrm{\Delta }\widehat{G}\right|},$$
(25)
where $`\lambda `$ is the constant defined through (14). The typical valley length is $`\widehat{L}1/H^{}`$, where $`H^{}`$ is given by (16). The valleys are spread out along the corridor at $`L\widehat{L}`$, separated by a typical distance $`M`$ which is also calculable from the large-deviation theory. For much larger $`L`$’s, the bubble energy becomes prohibitively high.

Clearly, the dynamics consists of two parts: At short times, it is dominated by the escape of the bubble out of an individual valley, and is analogous to the (biased) Sinai problem Sinai . At longer time scales, the bubble “hops” from one valley to another along the corridor of valleys. This dynamics, which is essentially that of a particle traversing a series of exponentially-distributed energy valleys \[see Eq. (25)\], has been extensively investigated previously in the context of the one-dimensional “trap” model slowa , slowb . Here we review some key results and refer the readers to Ref. berbou  for details. The basic quantity is the time $`\tau (\mathrm{\Delta }\widehat{G})e^{\beta \left|\mathrm{\Delta }\widehat{G}\right|}`$ to escape each valley of depth $`\mathrm{\Delta }\widehat{G}`$. The average time to traverse $`K`$ valleys over a length scale $`N=KM`$ by random walk is then given by

$$t_\mathrm{e}=K^2\tau N^2_{\lambda ^1}^{|\overline{E^{}}|}𝑑x\tau (x)𝒟(x),$$
(26)
where $`\tau `$ is the average of the trap time $`\tau (\mathrm{\Delta }\widehat{G})`$, and the limits of integration in (26) are from the magnitude of the typical valley depth $`\lambda ^1`$ to that of the deepest valley (13) expected for a segment of length $`N`$.

The total time according to Eq. (26) can be written as $`t_\mathrm{e}(N)N^z`$, with the dynamic exponent $`z`$ given by

$$z=\{\begin{array}{cc}2\hfill & \text{ for }\lambda >\beta \text{ (or }𝒯<𝒯_{\mathrm{loc}}\text{)}\hfill \\ 1+\beta /\lambda \hfill & \text{ for }\lambda <\beta \text{ (or }𝒯_{\mathrm{loc}}<𝒯<𝒯_\mathrm{d}\text{)}\hfill \end{array}$$
(27)
The anomalous exponent $`z>2`$ in the glass phase shows explicitly that the dynamics is slow, i.e., sub-diffusive.

### III.3 Glassy Dynamics

We next report the result of a Monte-Carlo simulation of the bubble dynamics on predefined random nucleotide sequences. We impose local dynamics in which the bubble can only change its length $`L`$ or shift its end position $`m`$ by a single base, as long as $`L1`$. To remove edge effects and probe the asymptotic dynamics, we use a very large sequence length ($`>10^4\mathrm{𝚋𝚙}`$) so that the bubble never reaches the boundary of the sequence given the duration of our numerical study. All disorder-averaged quantities reported are performed over $`10^4`$ random sequences.

1. Anomalous diffusion: To characterize the slow dynamics quantitatively, we show in Fig. 1(a) the time evolution of the average displacement $`R(t)=|m(t)m(0)|`$ of the bubble position for a few selected values of $`𝒯`$’s in the glass phase. Clearly, the displacement can be described by a power law of the form $`R(t)t^\nu `$, where we expect $`\nu =1/z`$. In Fig. 1(b), we plot the extracted exponents (circles) for different values of $`𝒯`$’s in the range $`𝒯_{\mathrm{loc}}𝒯<𝒯_\mathrm{d}`$. The expected values $`1/z`$ according to Eq. (27) (using the linear expression in (21) for $`\lambda `$) is shown as the solid line for comparison. We note that the observed exponents follow the general trend predicted, changing continuously from $`1/z=0.5`$ close to the expected location of the glass transition ($`𝒯_{\mathrm{loc}}0.8𝒯_\mathrm{d}`$), towards zero as $`𝒯𝒯_\mathrm{d}`$. For $`𝒯`$ close to $`𝒯_\mathrm{d}`$, the dynamics becomes exceedingly slow, making it difficult to access the asymptotic region. For $`𝒯𝒯_{\mathrm{loc}}`$, we also observed some finite-size effect. The overall agreement between the scaling theory and numerical results is within $`510\%`$ over the range tested.

In Fig. 2(a), we show the dependence of the average bubble length on time for different $`𝒯`$’s. The data depict the slow, logarithmic growth of the bubble length. Logarithmic growth is one of the signatures of glassy dynamics. Its occurrence in this particular system can be understood quantitatively as follows: The optimal bubble size $`L^{}(N)`$ in a segment of length $`N`$ depends logarithmically on $`N`$; see Eq. (15). On the other hand, for a bubble placed at an arbitrary position in a long sequence, the effective sequence length is the distance the bubble can explore within a time $`t`$, i.e., $`Nt^{1/z}`$ for the sub-diffusive dynamics expected in the glassy regime. Hence,

$$L^{}(t)\frac{1}{zH^{}}\mathrm{ln}t+\mathrm{const}.$$
(28)
is the expected length of the optimal bubble within a time t. Generally, we expect $`L^{}(t)`$ to be the upper bound of the observed bubble length $`L(t)`$, with $`L(t)L^{}(t)`$ for large $`t`$ deep in the glass phase. However below the glass transition, $`L(t)`$ must be finite even for $`t\mathrm{}`$.

In Fig. 2(b), we show the coefficients of the observed logarithmic time dependence of $`L(t)`$ for $`𝒯`$’s throughout the range $`𝒯_{\mathrm{loc}}<𝒯<𝒯_\mathrm{d}`$. Also shown is the upper bound $`1/(zH^{})`$ (solid line) according to (28), using the expression (22) for $`H^{}`$. We note that the difference between the data (circles) and the upper bound is nearly constant ($`1`$) for the range studied.

2. Aging: Perhaps the most characteristic feature of glassy dynamics is that the system “ages”, e.g., the temporal fluctuation of the system depends on how long the system has evolved from some (arbitrary) initial condition aging1 , aging : The longer it has evolved, the slower it fluctuates. This is easy to understand in the context of a rough energy landscape with deep valleys and high barriers, since the longer the system evolves, the deeper the energy valley it finds, and hence the higher the barrier it will have to overcome to travel farther. This feature is in marked contrast to sub-diffusive hydrodynamic systems which are time-translationally invariant.

Quantitatively, we can define the aging phenomenon via the time-dependent correlation function $`C(t_\mathrm{w},\mathrm{\Delta }t)`$, which measures how much the system changes in time $`\mathrm{\Delta }t`$, after first evolving for a waiting period $`t_\mathrm{w}`$ from the initial condition. Let us define a binary variable $`\eta _i(t)\{0,1\}`$, for each base $`i`$ of the nucleotide sequence. $`\eta _i(t)`$ takes on the value $`1`$ if base $`i`$ is open and belongs to the bubble at time $`t`$, and the value $`0`$ if base $`i`$ is paired. The correlation function, defined as $`C(t_\mathrm{w},\mathrm{\Delta }t)_i\eta _i(t_\mathrm{w})\eta _i(t_\mathrm{w}+\mathrm{\Delta }t)`$ after averaging over $`10000`$ random sequences, is a measure of the average self-overlap of the bubble at time $`t_\mathrm{w}`$ and $`t_\mathrm{w}+\mathrm{\Delta }t`$. A more convenient quantity to characterize is the fraction of overlap, $`C(t_\mathrm{w},\mathrm{\Delta }t)/L(t_\mathrm{w})`$, where $`L(t)=_i\eta _i(t)`$ is the instantaneous bubble length.

In Fig. 3(a), we show the overlap fraction, parameterized by the different waiting time $`t_\mathrm{w}`$’s for the system biased deep in the glass phase with $`𝒯=0.9𝒯_\mathrm{d}`$. The overlap fraction clearly depends on the waiting time, illustrating the glassy nature of the dynamics. In contrast, the same quantity computed for $`𝒯<𝒯_\mathrm{d}`$ (data not shown) gives no statistically significant dependence on $`t_\mathrm{w}`$. To characterize more quantitatively the behavior, we re-plot in Fig. 3(b) the curves in (a) with $`\mathrm{\Delta }t`$ normalized by $`t_\mathrm{w}`$. We find these curves to collapse reasonably onto a single master curve which exhibits a weak kink at $`\mathrm{\Delta }t/t_\mathrm{w}1`$. A naive explanation of this behavior is that for $`\mathrm{\Delta }tt_\mathrm{w}`$, the bubble stays approximately within the energy valley found at time $`t_\mathrm{w}`$, while for $`\mathrm{\Delta }tt_\mathrm{w}`$, the bubble makes excursion far away from the valley. For the one-dimensional “trap” model, it was shown rigorously corr-scale  that $`C(t_\mathrm{w},\mathrm{\Delta }t)`$ indeed scales as a function of $`\mathrm{\Delta }t/t_\mathrm{w}`$, even though the largest trap time actually scales sub-linearly with $`t_\mathrm{w}`$. This behavior can be understood in terms of the particle making multiple returns to the original valley after escaping it berbou , as manifested by the slow decay shown in Fig. 3(b) for $`\mathrm{\Delta }tt_\mathrm{w}`$.

## IV Discussion

In this study we investigated the thermodynamic and dynamic behaviors of twist-induced denaturation bubbles in a long, random sequence of DNA. The small bubbles associated with weak twist are delocalized, e.g., they flicker in and out of existence according to the Boltzmann distribution and are independent of the DNA sequence. The bubbles increase in lengths upon increase in the applied torque. When the largest bubbles reach a critical size $`L_{\mathrm{loc}}`$ which is of the order of a few tens of bases, the bubbles become localized to AT-rich segments which occur statistically in a long random sequence. According to the parameters SantaLucia  taken at $`37^{}𝙲`$ with $`[\mathrm{Na}^+]=1𝙼`$, the localization “transition” occurs at $`𝒯_{\mathrm{loc}}8\mathrm{pN}\mathrm{nm}`$, which is $`80\%`$ of the torque needed for bulk denaturation $`𝒯_\mathrm{d}`$. In the localized regime, the bubbles exhibit “aging” and move along the double helix sub-diffusively, with continuously varying dynamic exponents.

All of the results are obtained under the single-bubble approximation. Thermodynamically, we expect this approximation to be valid for DNA sequences of several thousand bases or less. This is due to the strongly cooperative nature of bubble formation, as manifested in the large initiation energy $`\gamma _1`$. The single bubble description of dynamics is further restricted by the finite life time of the bubble: Even at length scales where the single-bubble approximation is appropriate thermodynamically, the bubble may annihilate and reappear elsewhere in the sequence, effectively performing long-distance hops. Experimental knowledge of the bubble life time in the presence of an applied twist is needed to estimate the crossover time to the long-distance hopping regime. Qualitatively, we expect these bubbles to have much longer life times than the thermally denatured bubbles, since the applied twist plays the role of an energy barrier preventing bubble annihilation.

Finally, we note that bubble localization characterized in this study is a reflection of the statistical background present in long random nucleotide sequences. This background traps the bubble kinetically if the bubble size becomes sufficiently large. Thus, to localize denaturation bubbles at appropriate locations specified by designed sequences (e.g., promoters or replication origins) for biological functions, it is necessary to operate away from the localized regime, i.e., below the onset of localization.

This collaboration was made possible by the program on “Statistical Physics and Biological Information” hosted by the Institute for Theoretical Physics in Santa Barbara. The authors benefitted from discussions with D. Bensimon, R. Bundschuh, H. Chate, U. Gerland, D. Lubensky, M. Mezard, and Y.-k. Yu. TH is supported by NSF Grant No. 0211308 and a Burroughs-Wellcome functional genomics award. LT acknowledges the hospitality of UC San Diego where part of this work was carried out.


text:
# 1 Introduction

## 1 Introduction

One of the most important goals in the scientific program of a future linear electron–positron collider (LC) is to measure the properties of the Higgs boson precisely. Among them, the parity plays a prominent role. Depending on the mass and the Higgs mechanism, the Higgs parity can be measured by different means . In some cases, it can be best determined by investigating the properties of the Higgs boson decay into $`\tau `$ leptons . In  we studied the possibility of distinguishing a scalar from a pseudoscalar coupling of a light boson to fermions using its decay to a pair of $`\tau `$ leptons and their subsequent decays $`\tau ^\pm \rho ^\pm \overline{\nu }_\tau (\nu _\tau )`$ and $`\rho ^\pm \pi ^\pm \pi ^0`$. This decay chain gives an independent test of model and production mechanism. In fact, we have found that the angular distributions of the $`\tau ^\pm `$ decay products, which clearly distinguish the different parity states, are measurable using typical properties of a future detector at an $`e^+e^{}`$ linear collider. A special Monte Carlo technique was developed for that purpose . In the present paper we will extend that analysis, using the information that can be obtained from the measurement of the $`\tau `$ lepton impact parameter.

It turns out that the kinematic configuration of $`e^+e^{}\tau ^+\tau ^{}`$ events cannot be fully determined since the momenta of both $`\nu _\tau `$ and $`\overline{\nu }_\tau `$ are unknown. The lost information can be recovered partially by measuring energy and directions of all visible particles . When both $`\tau `$ leptons decay hadronically and all hadron momenta are recovered, the original $`\tau `$ direction can be obtained up to a twofold ambiguity . In the method, the requirement that the two $`\tau `$ leptons be essentially of equal energy and back to back in some known frame, is decisive. The measurement of the tracks of the hadrons, in particular their relative impact parameters, allows us to resolve this twofold ambiguity, see for instance .

In the case of linear colliders such as TESLA , a sufficiently precise reconstruction of the frame, where $`\tau `$ leptons are expected to be back to back turns out to be impossible, mainly due to beamstrahlung effect. This is a strong limitation for the Higgs boson parity measurement as proposed in Refs. . A way around was found in Ref. . Here, we will investigate if the measurement of the $`\tau `$ lepton impact parameter direction, with precision as expected for TESLA, can be used to optimize that method.

This work is organized as follows: in section 2 the geometry of $`\tau `$ production and decay is described, as well as our way to use the $`\tau `$ impact parameter in the measurement of the Higgs boson parity. In section 3 some details of the TAUOLA Monte Carlo simulation are given. Section 4 is devoted to the definitions used in our simulation of detector properties. In section 5 definitions of our observable and numerical results are given. The summary, section 6, closes the paper.

## 2 Description of the method

In the decay of the $`\tau `$ lepton into one charged particle, the $`\tau `$ impact parameter can be defined as the distance of closest approach (in the plane perpendicular to the beam axis) of the charged particle to a reference point, which is assumed to be the production point of the $`\tau `$ (see Fig. 1). A positive sign is given to this quantity if the track crosses the $`\tau `$ direction in the same hemisphere, relative to the $`\tau `$ production point, in which it lays, and negative otherwise. The production point is assumed to lie on the beam line. In our study we will exploit the possibility to determine the direction of the impact parameter, rather than the related distance. In fact we will be interested, where the $`\tau `$ production point is localized with respect to the charged pion track.

Before we address the question of the information available from the impact parameter, let us recall some details of the actual observable proposed in Ref. for the Higgs boson parity measurement. The method relies on measuring the acoplanarity angle of the two planes, spanned on $`\rho ^\pm `$ decay products and defined in the $`\rho ^+\rho ^{}`$ pair rest frame; some additional selection cuts also need to be applied. Let us present the elements of the observable as the following points.

1. Acoplanarity of the $`\rho ^\mathbf{+}`$ and $`\rho ^{\mathbf{}}`$ decay planes:

* The reconstructed four-momenta of $`\pi ^+`$ and $`\pi ^0`$ are combined to yield the $`\rho ^+`$ four-momentum. The same is done for the $`\rho ^{}`$.
    * All reconstructed four-momenta are boosted into the $`\rho ^+\rho ^{}`$ rest frame.
    * The angle $`\phi `$ between the planes of the $`\rho ^+`$ and $`\rho ^{}`$ decay products in this frame is the acoplanarity.
2. Normalized energy differences:
The events are divided into two classes depending on the value of $`y_1y_2`$, where

$$y_1=\frac{E_{\pi ^+}E_{\pi ^0}}{E_{\pi ^+}+E_{\pi ^0}};y_2=\frac{E_\pi ^{}E_{\pi ^0}}{E_\pi ^{}+E_{\pi ^0}}.$$
(1)
The energies of $`\pi ^\pm ,\pi ^0`$ are to be taken in the respective $`\tau ^\pm `$ rest frames.
3. Replacement $`\tau `$ rest frame:
In Ref.  the method of reconstruction of replacement $`\tau ^\pm `$ rest frames was proposed. We will use this method here as well, but we will omit the details. In fact they are rather unimportant. The energies of pions taken directly from the laboratory frame can equally well be used in the above definitions of Eq. (1).
4. Higgs rest frame:
If the information on the beam energies and the energies of all other observed particles is taken into consideration, the Higgs rest frame can be reconstructed, for instance in the Higgsstrahlung production process, $`e^+e^{}HZ`$, when the $`Z`$ boson decays either into a charged lepton pair or hadronically. We may define the ‘reconstructed’ Higgs boson momentum as the difference of the sum of beam momenta and momenta of all visible particles, i.e. decay products of $`Z`$ and all radiative photons of $`|\mathrm{cos}\theta |<0.98`$.
5. $`𝝉`$ energy:

* In the reconstructed Higgs boson rest frame, the $`\tau `$ four-momenta are estimated in a crude way by assuming the direction of the respective $`\rho `$’s and an energy of $`m_H/2`$.
    * The $`\tau `$ momenta are boosted back to the laboratory frame and their energies are taken to be the reconstructed energies of $`\tau `$ leptons.
6. $`𝝉`$ direction:
The $`\tau `$ direction in the laboratory system is constrained by two requirements:

* It has to lie on a cone around the $`\rho `$ direction with opening angle $`\psi `$ defined as:

    $$\mathrm{cos}\psi =\frac{2E_\tau E_\rho m_\tau ^2m_\rho ^2}{2\stackrel{}{p}_\tau \stackrel{}{p}_\rho },$$
    (2)
    where $`E`$, $`\stackrel{}{p}`$, $`m`$ are the energy, momentum and mass, respectively. The four-momentum of $`\rho `$ is measured as the sum of the four-momenta of its charged and neutral daughters.
    * It has to lie in the plane spanned by the vector pointing from $`PP`$ to $`PCA`$ and the $`\pi ^+`$ momentum.

The intersection of the cone and the plane is calculated numerically. There are three cases:

* One solution: the solution is taken as the $`\tau `$ direction (see Fig. 2 case A).
    * Two solutions: one of the two solutions (see Fig. 2 case B) is taken on a random basis.
    * No solution: the direction on the cone closest to the $`PCA`$$`\pi ^+`$ momentum plane is taken (see Fig. 2 case C).
7. Impact parameter improved replacement $`𝛕`$-rest frame:
With the help of the $`\tau `$ energy and $`\tau `$ momentum defined in points 5 and 6 above, a new impact parameter improved replacement $`\tau `$ rest frame can be defined. In this way we have an alternative way of estimating the difference of $`\pi ^\pm `$, $`\pi ^0`$ energies see Eq. (1) in $`\tau ^\pm `$ rest frames.

## 3 Monte Carlo set-up

In the following discussion all the Monte Carlo samples have been generated with the TAUOLA library . For the production of $`\tau `$ lepton pairs, the Monte Carlo program PYTHIA 6.1  was used111 It was shown that the interface can work in the same manner with the HERWIG  generator.. The production process $`e^+e^{}ZH\mu ^+\mu ^{}(q\overline{q})H`$ has been chosen, with a Higgs boson mass of $`120`$ GeV and a centre-of-mass energy of $`350`$ GeV.

The effects of initial-state bremsstrahlung were included in the PYTHIA generation. For the $`\tau ^\pm `$ lepton pair decay (with full spin effects from the $`H\tau ^+\tau ^{}`$, $`\tau ^\pm \rho ^\pm \overline{\nu }_\tau (\nu _\tau )`$, $`\rho ^\pm \pi ^\pm \pi ^0`$ chain) the interface explained in Ref.  was used. It is an extended version of the interface of Refs. .

## 4 Detection parameters

To test the feasibility of the measurement, some assumptions about the detector effects have to be made. We include, as the most critical for our discussion, effects due to inaccuracies in the measurements of the $`\pi ^\pm ,\pi ^0`$ momenta and $`\tau ^\pm `$ leptons impact parameters. We assume Gaussian spreads of the ‘measured’ quantities with respect to the generated ones, and we use the following algorithm to reconstruct the energies of $`\pi `$’s in their respective $`\tau ^\pm `$ rest frames.

1. Charged-pion momentum: We assume a 0.1% spread on its energy and direction.
2. Neutral-pion momentum: We assume an energy spread of $`\frac{5\%}{\sqrt{E[\mathrm{GeV}]}}`$. For the $`\theta `$ and $`\varphi `$ angular spread we assume $`\frac{1}{3}\frac{2\pi }{1800}`$. (These neutral pion resolutions can be achieved with a 15% energy error and a 2$`\pi `$/1800 direction error in the gammas resulting from the $`\pi ^0`$ decays.) These resolutions have been approximately verified with a SIMDET , a parametric Monte Carlo program for a TESLA detector .
3. The reconstructed Higgs boson rest frame: We assume a spread of $`2`$ GeV with respect to the transverse momentum of the Higgs boson, and $`5`$ GeV (because of the beamstrahlung effect) for the longitudinal component.
4. The impact parameter:
The angular resolution of the vector pointing from $`PP`$ to $`PCA`$ (see Fig. 1) has been simulated for the TESLA-like detector. The simulation is based on the anticipated performance of a 5-layer CCD vertex detector as described in . For Higgsstrahlung events with $`H\tau ^+\tau ^{}`$ and $`\tau ^\pm \rho ^\pm \overline{\nu }_\tau (\nu _\tau )`$ at $`m_H=`$ $`120`$ GeV and $`\sqrt{s}=`$ $`350`$ GeV, the angular resolution has been found to be approximately 25. We use this resolution in our simulations.

## 5 Acoplanarity of the $`\rho ^+`$ and $`\rho ^{}`$ decay products with the help of the impact parameter

For all figures in our paper we will use the same angle $`\phi `$ built from smeared four-momenta of the $`\rho ^\pm `$ decay products. The usefulness of this acoplanarity distribution manifests itself only after selection cuts are applied.

In Fig. 3, our $`\rho ^+\rho ^{}`$ decay products acoplanarity distribution angle $`\phi `$, defined in the rest frame of the reconstructed $`\rho ^+\rho ^{}`$ pair is shown. Unobservable generator-level $`\tau ^\pm `$ rest frames are used for the calculation of selection cuts (Eq. (1)). The two plots represent events selected by the differences of $`\pi ^\pm `$ $`\pi ^0`$ energies defined in these respective $`\tau ^\pm `$ rest frames. On the left plot, the signs are required to be the same $`y_1y_2>0`$, whereas on the right one, events are taken with $`y_1y_2<0`$. This figure quantifies the size of the parity effect in an idealized condition, which we will attempt to approach with realistic ones. The size of the effect was substantially diminished when the detector set-up was included for $`\tau ^\pm `$ rest frames reconstruction as well, see Fig. 4, following exactly the method of Ref. . Let us use this result later, as a reference point to quantify the size of our improvements.

It turns out that a slight improvement is visible, when the impact parameter is taken for the $`\tau ^\pm `$ rest frames determination, as explained in section 4. In fact, the maximum relative difference between the scalar and pseudoscalar cases is larger by $`12\%`$ (compare Figs. 4 and 5) if a Gaussian spread of the $`\tau `$ impact parameter $`\sigma _{IP}=25^{}`$ is used. We have also checked that the precision of the $`\tau `$ lepton impact parameter measurement is not critical in this case. The twofold ambiguities as explained in Fig. 2 are the main reason for the sensitivity to be much lower than in Fig. 3.

At the cost of introducing cuts and thus reducing the number of accepted events, we can go one step further in the improvement of the method. We first start by requiring that the signs of the energy differences $`y_1`$ and $`y_2`$ be the same with both methods (with and without $`\tau `$ lepton impact parameter). The result is depicted in Fig. 6. Only $`52\%`$ of events are accepted, but the parity effect grows by $`107\%`$ in comparison to Fig. 4.

The events can also be treated differently, according to the number of solutions; see Fig. 2. If the same cut is imposed only for $`\tau ^\pm `$ leptons with double solutions, while for single-solution events the $`\tau ^\pm `$ rest frames obtained with the help of the impact parameter are used alone, the number of accepted events turns out to be $`72\%`$, but the parity effect is larger by only $`67\%`$ than in Fig. 4. This result is depicted in Fig. 7.

Other cuts are also envisageable, but the two proposed here are easy to implement and lead to an acceptable loss of statistics, together with a sizeable increase of parity effect so that there is not much room for further improvements.

## 6 Summary

We have found that the measurement of the impact parameter can be helpful in the determination of the Higgs boson parity at a future linear collider, using the $`H/A\tau ^+\tau ^{}`$; $`\tau ^\pm \rho ^\pm \overline{\nu }_\tau (\nu _\tau )`$ cascade decay. Thanks to the improved method the effect of parity is enhanced by $`107\%`$ at the cost of losing only $`48\%`$ events, which represents a gain of a factor of $`1.5`$ in sensitivity with respect to the analysis without impact parameter. This correspond to $``$ 4.5$`\sigma `$ separation between scalar and pseudoscalar for the $`120`$ GeV Standard Model Higgs and other assumptions similar to the one taken in Ref. .

We have also checked, that reasonable variations in the assumptions of detection uncertainties in measurement of charged and neutral $`\pi `$’s as well as of the $`\tau ^+`$ and $`\tau ^{}`$ impact parameters, hardly change the sensitivity of our method at all. This is because the dominant uncertainty is saturated by the twofold ambiguity as explained in Fig. 2. We are afraid that a reconstruction of the Higgs boson momentum to precision better than the $`\tau `$ mass would be needed to resolve the ambiguity. This seems unrealistic for several reasons, e.g. the beamstrahlung effect.

Probably, a fit to the multidimensional distribution might be more efficient than the simple cuts applied in this work. We believe, however, that such studies should be done only when the detector properties have been understood better. Finally, we expect that further improvements, e.g. involving the $`\tau ^\pm a_1^\pm \nu _\tau `$ decay mode, may turn out to be helpful for the measurement of the Higgs boson parity.

### Acknowledgements

This work is partly supported by the Polish State Committee for Scientific Research (KBN) grants Nos 5P03B09320, 2P03B00122 as well as the BMBF (WTZ) project number POL 01/103 and by the DESY $`\&`$ INP Cracow collaboration grant TESLA.

One of the authors (M.W.) is grateful for the warm hospitality extended to her by the FLC group at DESY in Hamburg. She also would like to thank the “Marie Curie Programme” of the European Commission for a fellowship.

Useful discussions with J. H. Kühn are also acknowledged.