text:
# Is there a d.c. Josephson effect in bilayer quantum Hall systems?

## Abstract

We argue on the basis of phenomenological and microscopic considerations that there is no d.c. Josephson effect in ordered bilayer quantum Hall systems, even at $`T=0`$. Instead the tunnel conductance is strongly enhanced, approaching a finite value proportional to the square of the order parameter as the interlayer tunneling amplitude vanishes.

Introduction: Weakly disordered bilayer quantum Hall systems with small inter-layer separation $`d`$ have broken symmetry ground states that can be regarded either as as easy-plane ferromagnets  or as excitonic superfluids. Experimental evidence for this state, in which phase coherence between different layers occurs spontaneously, was first uncovered in quantum Hall effect activation energy studies  that showed anomalous in-plane field dependence indirectly  related to its collective properties. Recently, however, Spielman et al.  discovered a spectacular feature in the zero-bias conductance $`G_T`$ that appears when bilayer quantum Hall system parameters are tuned into the regime where order occurs. Although it appears clear that the dramatic conductance peak they observe is due to collective inter-layer tunneling in the ordered state, it has not yet been possible to explain its height and width or its dependence on proximity to the ordered state phase boundary. In particular, the close similarity between the phenomenological effective theory  of an ordered bilayer quantum Hall system and that of a Josephson junction suggests  that bilayers should exhibit a d.c. Josephson effect, i. e., that persistent currents can flow between the layers without any bias potential. Recent theoretical work  has addressed the tunneling characteristics of clean and disordered bilayer quantum Hall systems at finite bias voltages, finite in-plane fields, and finite temperatures, demonstrating among other things that in the ordered state the inter-layer tunneling amplitude $`\mathrm{\Delta }_t`$ cannot be treated perturbatively, i.e. that $`G_T/\mathrm{\Delta }_t^2`$ diverges. The central question concerning these experiments, whether or not a d.c. Josephson effect occurs in principle, has not yet been addressed directly, although divergent views have been expressed by different authors. Experimentally there is no evidence for a d.c. Josephson effect, i.e. the $`T0`$ zero-bias conductance peak appears to be finite. The explanation for this finding need not be fundamental, however; for example the $`T0`$ order-parameter could vanish  in current samples due to quantum fluctuations that are expected to be enhanced by disorder, or the highest measured conductance values could be limited by extrinsic experimental effects. In this Letter we argue that in bilayer quantum Hall systems $`G_T`$ always remains finite at $`T=0`$. Nevertheless, $`G_T/\mathrm{\Delta }_t^2`$ does diverge for $`\mathrm{\Delta }_t0`$ and inter-layer tunneling cannot be treated perturbatively. In the following paragraphs, we first discuss the physical picture that underlies our theory and then present an approximate but fully microscopic calculation that we believe captures all essentials of the effect.

Phenomenological Theory: We use a pseudospin language to describe bilayer quantum Hall systems.

$$S_\alpha =\frac{1}{2}\underset{k,\sigma ^{},\sigma }{}c_{k\sigma ^{}}^{}\tau _{\sigma ^{}\sigma }^{(\alpha )}c_{k\sigma }$$
(1)
is the total pseudospin component in direction $`\alpha `$, $`k`$ is a Landau level orbital labels, $`\sigma ,\sigma ^{}`$ are the pseudospin labels, and $`\tau _{\sigma ^{}\sigma }^{(\alpha )}`$ are Pauli spin matrices with pseudospin up/down representing electrons in top/bottom layers. In this language the ordered state is an $`XY`$ easy-plane pseudospin ferromagnet. Its long-wavelength, low-frequency, small-amplitude dynamics should be described  by the linearized Landau-Lifshitz-Gilbert equation  of magnetism

$`{\displaystyle \frac{dM_z}{dt}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }_t}{\mathrm{}}}M_y,`$ (2)
$`{\displaystyle \frac{dM_y}{dt}}`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{}}}\left[eV(\mathrm{\Delta }_t+4\pi l^2\beta )M_z\right]{\displaystyle \frac{M_y}{\tau }}.`$ (3)
In Eq.( 3) $`\widehat{M}(1[M_y^2+M_z^2]/2,M_y,M_z)`$ is a unit vector which specifies the pseudospin ordered moment orientation, and $`\beta `$ is the pseudospin anisotropy energy per unit area  which will be discussed at greater length below. For quantum Hall ferromagnets there is no relaxation term in the first of Eqs.( 3) because only inter-layer tunneling, which we take to be constant, violates the separate conservation of charge in each layer. The factor in square brackets in the second of Eqs.( 3) is the effective field in the $`\widehat{z}`$-direction which includes both external and induced contributions to the electrochemical potential difference $`\delta \mu `$ between the two layers. It follows from the first of the Eqs.( 3) that the inter-layer current is $`I=(NeM_0/2)(dM_z/dt)=(NeM_0\mathrm{\Delta }_t/2\mathrm{})M_y`$, where $`N`$ is the total number of electrons, $`lim_{\mathrm{\Delta }_t0}S_x=NM_0/2`$ is the pseudospin ordered moment. Here $`M_0`$ is a dimensionless order parameter that approaches $`1`$ for layer separation $`d0`$ in the absence of disorder . In the steady state, the driving and relaxation terms in the second of Eqs.( 3) cancel, from which it follows that $`M_y=\delta \mu \tau /\mathrm{}`$ and that the tunnel conductance is

$$G_T=\frac{e^2}{\mathrm{}}\frac{NM_0\mathrm{\Delta }_t\tau }{2\mathrm{}}.$$
(4)
The relaxation time in bilayer quantum Hall systems can be evaluated microscopically by evaluating the dynamic pseudospin response function. It follows from Eq.( 3) that at low frequencies

$$\chi _{yz}(\omega )=\frac{M_0i\mathrm{}\omega }{\left[(\mathrm{\Delta }_t+4\pi l^2\beta )\mathrm{\Delta }_t\mathrm{}^2\omega ^2\right]i\mathrm{}^2\omega /\tau }.$$
(5)
We derive a response function of precisely this form from a microscopic theory below.

Is there a d.c. Josephson Effect? In using Eq.( 4) we are asserting that the d.c. conductance is equal to the $`\omega 0`$ limit of the a.c. conductance and thus denying the possibility of persistent d.c. tunnel currents. In a conventional Josephson junction geometry, persistent currents are enabled by order parameter phase rigidity across the entire system. From the point of view of microscopic mean-field-theory, the anomalous self-energy couples electrons and holes, forming an equilibrated quasiparticle system which self-consistently establishes different order parameter phase values on opposite sides of the junction. In other words, the Dyson equation has a set of self-consistent solutions with continuously variably order parameter phase difference across the junction. In a quantum Hall system, on the other hand, there is no analog of overall phase rigidity, only phase difference rigidity supported by non-local interlayer interactions. From the microscopic point of view, the exchange self-energy associated with order couples electrons in balanced layers through a pseudospin effective field that points in the same direction as the pseudospin order parameter. Since inter-layer tunneling adds a pseudospin effective field in the $`\widehat{x}`$ direction, a self-consistent equilibrated quasiparticle system is possible only when the pseudospin order points in the $`\widehat{x}`$ direction, in which case no current flows. A collective interlayer current can always decay by making particle-hole transitions in the non-equilibrium quasiparticle system. The microscopic calculations presented below describe this effect. Of course, persistent currents can also decay in standard Josephson junctions, even at $`T=0`$, because of collective quantum tunneling of the order parameter field. This decay mechanism is, however, qualitatively weaker and easily distinguished.

Microscopic Theory: To describe a quantum-Hall bilayer, we use a self-consistent Born approximation for the disorder self-energy, a self-consistent Hartree-Fock approximation for the interaction self-energy, and include consistent ladder and bubble vertices in the two-particle Greens function of interest. It is essential for our analysis that this approximation, summarized diagrammatically in Fig. 1, captures the physics associated with the separate conservation of charge in each layer when $`\mathrm{\Delta }_t0`$ and that all necessary diagram sums can be evaluated accurately. The disorder averaged Greens function depends on frequency only and satisfies the Dyson equation

$$𝒢^1(i\nu _n)=\left(\begin{array}{cc}i\nu _n+\rho _{TT}\mathrm{\Gamma }_A+v_s^2𝒢_{TT}(i\nu _n)& +\mathrm{\Delta }_t/2+\rho _{TB}\mathrm{\Gamma }_E+v_d^2𝒢_{TB}(i\nu _n)\\ +\mathrm{\Delta }_t/2+\rho _{BT}\mathrm{\Gamma }_E+v_d^2𝒢_{BT}(i\nu _n)& i\nu _n+\rho _{BB}\mathrm{\Gamma }_A+v_s^2𝒢_{BB}(i\nu _n)\end{array}\right),$$
(6)
where $`\mathrm{\Gamma }_A`$ and $`\mathrm{\Gamma }_E`$ are intralayer and interlayer exchange  integrals , $`v_s^2`$ and $`v_d^2`$ are the disorder correlation functions in same and different layers, and the labels $`T`$ and $`B`$ refer to top and bottom layers. In Eq.( 6) $`\rho _{\sigma \sigma ^{}}`$ is the density matrix which is determined self-consistently by integrating Greens function spectral weights up to the Fermi energy. We concentrate here on balanced bilayer systems at total filling factor $`\nu =1`$ so that $`\rho _{TT}=\rho _{BB}=1/2`$, and take $`v_d^2=0`$ since the disorder potentials in different layers are not expected to be correlated. $`\rho _{TB}=\rho _{BT}^{}`$ is non-zero when order is established. Unlike the analogous Josephson-junction-system Dyson equation, when $`\mathrm{\Delta }_t0`$, Eq.( 6) has a solution only for purely real $`\rho _{TB}`$. Persistent currents would occur if complex solutions existed.

Inversion symmetry in balanced bilayers separates the $`y`$ and $`z`$ pseudospin response function components, related to correlations between many-particle states with opposite parity, from the charge and $`x`$ pseudospin response functions, related to correlations between states with the same parity. A somewhat lengthy but elementary calculation that follows line similar to earlier work  leads to the following expressions:

$$\left[\begin{array}{cc}\chi _{yy}(\omega )& \chi _{yz}(\omega )\\ \chi _{zy}(\omega )& \chi _{zz}(\omega )\end{array}\right]^1=\left[\begin{array}{cc}\mathrm{\Gamma }_E& 0\\ 0& 2V_x\mathrm{\Gamma }_A\end{array}\right]+\left[\begin{array}{cc}\mathrm{\Pi }_{yy}(\omega )& \mathrm{\Pi }_{yz}(\omega )\\ \mathrm{\Pi }_{zy}(\omega )& \mathrm{\Pi }_{zz}(\omega )\end{array}\right]^1,$$
(7)
where $`2V_x`$ is the difference between intra-layer and inter-layer Coulomb interactions,

$`\mathrm{\Pi }_{\alpha \beta }(i\omega )`$ $`=`$ $`{\displaystyle \frac{1}{\beta }}{\displaystyle \underset{i\nu _n}{}}\left[𝒮^1v\right]_{\alpha \beta }^1(i\omega ,i\nu _n),`$ (8)
$`𝒮_{\alpha \beta }`$ $`=`$ $`{\displaystyle \underset{\sigma _1\sigma _2\sigma _1^{^{}}\sigma _2^{^{}}}{}}\tau _{\sigma _1^{^{}}\sigma _1}^{(\alpha )}𝒢_{\sigma _2\sigma _1^{^{}}}(i\nu _n)𝒢_{\sigma _1\sigma _2^{^{}}}(i\omega +i\nu _n)\tau _{\sigma _2^{^{}}\sigma _2}^{(\beta )},`$ (9)
and the only non-zero matrix element of $`v`$ is $`v_{zz}=v_s^2`$.

In this generalized RPA theory, the $`\mathrm{\Pi }`$ response functions are those of non-interacting electrons with uncorrelated random potentials in the two layers and an effective tunnel splitting $`\mathrm{\Delta }_{qp}=\mathrm{\Delta }_t+2\rho _{TB}\mathrm{\Gamma }_E`$ which we refer to as the quasiparticle tunneling amplitude. It follows from Eqs.( 8) and ( 9) that at low frequencies

$$\left[\begin{array}{cc}\mathrm{\Pi }_{yy}(\omega )& \mathrm{\Pi }_{yz}(\omega )\\ \mathrm{\Pi }_{zy}(\omega )& \mathrm{\Pi }_{zz}(\omega )\end{array}\right]^1=\left[\begin{array}{cc}\mathrm{\Delta }_{qp}/M_{qp}& +i\mathrm{}\omega /M_{qp}\\ i\mathrm{}\omega /M_{qp}& \mathrm{\Delta }_{qp}/M_{qp}i\mathrm{}^2\omega /\mathrm{\Delta }_{qp}\tau _{qp}M_{qp}\end{array}\right]$$
(10)
where the factor of $`\mathrm{}\omega `$ in the damping term of the $`zz`$ component of this matrix reflects the phase space available at $`T=0`$ for absorption of collective motion energy by the quasiparticle system. Solving for $`\mathrm{\Pi }_{yz}`$ we find

$$\mathrm{\Pi }_{yz}(\omega )=\frac{M_{qp}i\mathrm{}\omega }{\left[\mathrm{\Delta }_{qp}^2\mathrm{}^2\omega ^2\right]i\mathrm{}^2\omega /\tau _{qp}},$$
(11)
the non-interacting electron version of the phenomenological Landau-Lifshitz-Gilbert expression. Eq.( 5) and Eq.( 11) motivate the notation chosen for the coefficient of $`i\mathrm{}\omega `$ in $`\left(\mathrm{\Pi }^1\right)_{zz}`$ in Eq.( 10).

The tunnel conductance for this fictional non-interacting electron system is

$$G_{qp}=\frac{e^2}{\mathrm{}}\frac{NM_{qp}\mathrm{\Delta }_{qp}\tau _{qp}}{2\mathrm{}}.$$
(12)
For small $`\mathrm{\Delta }_{qp}`$, it follows from the standard Golden rule argument that $`G_T\mathrm{\Delta }_{qp}^2`$ which, since $`M_{qp}`$ is proportional to $`\mathrm{\Delta }_{qp}`$, implies that $`\tau _{qp}`$ approaches a constant as $`\mathrm{\Delta }_{qp}0`$. This result is expected since in this limit, the two quasiparticle systems are independent and $`\mathrm{\Pi }_{zz}(\omega )=(\mathrm{\Delta }_t/i\mathrm{}\omega )\mathrm{\Pi }_{yz}(\omega )`$ is the standard single-layer diffusive response response.

Combining Eq.( 7) and Eq.( 10) we obtain

$$\left[\begin{array}{cc}\chi _{yy}(\omega )& \chi _{yz}(\omega )\\ \chi _{zy}(\omega )& \chi _{zz}(\omega )\end{array}\right]^1=\left[\begin{array}{cc}\mathrm{\Delta }_t/M_{qp}& +i\mathrm{}\omega /M_{qp}\\ i\mathrm{}\omega /M_{qp}& \left[\mathrm{\Delta }_t+\beta \right]/M_{qp}i\mathrm{}^2\omega /\mathrm{\Delta }_{qp}\tau _{qp}M_{qp}\end{array}\right],$$
(13)
where we have used that $`\mathrm{\Delta }_{qp}=\mathrm{\Delta }_t+M_{qp}\mathrm{\Gamma }_E`$ and recalled  that $`\beta =M_{qp}(2V_x+\mathrm{\Gamma }_A\mathrm{\Gamma }_E)`$ is the GRPA theory result for the anisotropy energy per unit area. In Fig.( 2) we plot $`\mathrm{}\chi _{yz}^1(\omega )`$ vs. $`\omega `$ for various values of $`\mathrm{\Delta }_t,v_s^2`$, and compare the SCBA response functions with Eq.( 5). The inset to Fig.( 2) shows the dependence of the quasiparticle relaxation rate $`\tau _{qp}`$ on $`\mathrm{\Delta }_{qp}`$ obtained by fitting to Eq.( 11). It follows that GRPA theory recovers the phenomenological Eq.( 5) with $`M_0=M_{qp}=2\rho _{TB}`$ and relaxation time $`\tau =\tau _{qp}\mathrm{\Delta }_{qp}/\mathrm{\Delta }_t`$. The GRPA conductance

$$G_T=\frac{e^2}{\mathrm{}}\frac{NM_{qp}\mathrm{\Delta }_{qp}\tau _{qp}}{2\mathrm{}}=G_{qp}.$$
(14)
Discussion: Our microscopic theory predicts that in the ordered state $`G_T`$ remains finite as $`\mathrm{\Delta }_t0`$, consistent with the extremely large enhancements seen experimentally. At small $`\mathrm{\Delta }_t`$, the experimental situation, Eq.( 14) predicts that $`G_T\mathrm{\Delta }_{qp}M_{qp}M_{qp}^2`$, consistent with the smooth monotonic growth of the conductance peak seen experimentally after the ordered state is established. We emphasize that the absence of a Josephson effect in the bilayer tunneling conductance is not related to superflow within the 2D planes; these bilayer systems are 2D excitonic  superfluids and can carry oppositely directed currents in the two layers without dissipation. As pointed out by Wen and Zee  an analog of the Josephson effect for these counter propagating currents can be realized by creating weak links within the 2D layers. We note that the ordered state tunnel conductance in our theory can grow to values larger than the in-plane bilayer conductivity, suggesting that current-path and condensate to quasiparticle current conversion issues , as well as the spatial inhomogeneities known to exist on quantum Hall plateaus , could play a role in general in the interpretation of side-contacted tunneling studies .

This work was supported in part by the Welch Foundation and by the Indiana 21st Century Fund. The authors acknowledge helpful conversations with Anton Burkov, Jim Eisenstein, Steve Girvin, Leo Radzihovsky, John Schliemann, and Ady Stern.


format: text

text:
# Abelian threshold models and forced weakening

## Abstract

Mean field slider block models have provided an important entry point for understanding the behavior of discrete driven threshold systems. We present a method of constructing these models with an arbitrary frictional weakening function. This ‘forced weakening’ method unifies several existing approaches, and multiplies the range of possible weakening laws. Forced weakening also results in Abelian rupture propagation, so that an avalanche size depends only on the initial stress distribution. We demonstrate how this may be used to accurately predict the long-time event statistics of a simulation.

05.45.Ra, 07.05.Tp, 05.65.+b

Complex spatial or temporal patterns often emerge when systems are driven from equilibrium and experience instabilities. Coupled lattice maps have become a popular method of modelling such phenomena, seeking to capture the essential physics in a discrete form amenable to numerical simulation. In the study of earthquakes, simple slider-block models  have long been employed to investigate the origin of magnitude-frequency scaling. Similar models are used to describe such diverse phenomena as pinned charge density waves , flux lattices in type II superconductors , and creeping contact lines .

Recently, analytic results for the avalanche size distribution have been presented for mean-field slider-block models . These models are characterized by homogeneous, infinite range elastic interactions, and exhibit complex event histories and regimes of behavior. Despite their simplicity, mean-field models remain sensitive to the choice of update rules and details of implementation. This is especially evident in models that impose a particular form of simulated frictional weakening , where different modes of behavior appear as the strength of the weakening is varied.

Here we describe a method of constructing a mean-field threshold model where one may simulate arbitrary weakening laws under identical rules of evolution. This unifies the analysis of previously incompatible models, and provides more freedom in numerical simulation. This technique also results in Abelian rupture propagation, where the size of a simulated earthquake is uniquely determined from initial conditions, leading to a rigorous and implementation-independent analysis.

The general slider-block model represents stick-slip motion along a fault plane with $`N1`$ discrete coordinates (or ‘sites’). Each site $`i`$ is assigned a slip deficit $`u_i`$, which measures the local distance from elastic equilibrium. The sites are pinned in place by frictional forces, and are subject to a restoring force (stress) proportional to their slip deficit. Internal disorder gives rise to an additional component of stress due to elastic interactions. The stress $`\stackrel{~}{\sigma }_i`$ at a site $`i`$ is related to the slip deficits through a linear constitutive relation

$$\stackrel{~}{\sigma }_i=K_Lu_i\underset{j}{}K_{ij}\left(u_iu_j\right)$$
(1)
where $`K_L`$ and $`K_{ij}`$ are spring constants. If we impose uniform (mean-field) interactions between all the elements, $`K_{ij}=K_C/N`$, the above relation becomes

$$\stackrel{~}{\sigma }_i=K_Lu_iK_C(u_iu)$$
(2)
where $`u=N^1_iu_i`$ will denote an average over all $`N`$ sites in the model. We obtain a unitless expression by dividing by $`K_Ca`$, where $`a`$ is a characteristic microscopic length. Defining the unitless slip deficit $`\varphi =u/a`$, stress $`\sigma =\stackrel{~}{\sigma }/(K_ca)`$, and spring constant ratio $`K_R=K_L/K_C`$, Eq.  simplifies to

$$\sigma _i=(K_R+1)\varphi _i+\varphi .$$
(3)
For finite $`N`$ we will refer to this as the near mean-field (NMF) model. Note that it is easy to invert Eq. (3) for the slip deficits in terms of stresses,

$$\varphi _i=\frac{\sigma _i}{K_R+1}\frac{\sigma }{K_R(K_R+1)}$$
(4)
so the configuration is uniquely determined by the parameter $`K_R`$ and either the slip deficits or stresses alone.

The model is slowly driven away from equilibrium by uniformly increasing the slip deficits. Eventually the stress at one site will surpass the maximum local frictional force and ‘fail’, sliding toward its equilibrium point. The motion of a failed site will change the mean slip deficit $`\varphi `$, and produce a change in stress at other sites. If this change brings other sites to their threshold, they will also fail, producing an avalanche interpreted as a single event.

It is assumed that following the initiation of motion, the frictional force on a site will weaken, producing a transient dynamic instability. In discrete time we cannot model the dynamic slip or velocity of the site, but instead assign a residual stress $`\sigma ^R`$ at which the motion arrests. This $`\sigma ^R`$ is chosen from a probability distribution independently for each failed site. Since slips occur instantaneously, we lose the interplay between a continuously evolving stress field and frictional force at a site. The behavior of dynamical models is known to strongly depend on the form of frictional weakening , so the ability to include equivalent effects in discrete models would be advantageous.

Since the NMF model is not dynamical we are only interested in large-scale features of its behavior that are independent of microscopic dynamics. Thus we are free to choose the simplest update rules that are consistent with the phenomenon of interest. In practice, we assume that a single site $`j`$ reaches its stress threshold first. Since the physics will depend only on changes in stress, we may impose a uniform failure threshold $`\sigma ^F`$ by absorbing any variations into the residual stress distribution. The slip displacement $`\mathrm{\Delta }_j=\varphi _j^{(f)}\varphi _j^{(i)}`$ is related to the change in stress $`\mathrm{\Delta }\sigma _j=\sigma _j^R\sigma ^F`$ by $`\mathrm{\Delta }_j=\mathrm{\Delta }\sigma _j/(K_R+1N^1)`$. The motion of the site will change the mean slip deficit $`\varphi `$ by $`\mathrm{\Delta }_j/N`$. We may view this as a transfer of stress from failing sites to all others.

The above describes a series type dislocation where sites fail in sequence and the stress transfer occurs to other sites all at once. This makes it likely that any failing site (other than the single initiator) will have a stress slightly above the threshold, which subtly provides an order-of-failure dependence to the stress transfer. As a consequence, the exact stress transfer in simulation will depend on obscure factors like the order of iteration over sites. To eliminate this we must examine the stress transfer in more detail.

Suppose that in the course of a (possibly ongoing) event there have been $`k`$ block failures. Let $`\{k\}`$ represent the set of indices of failed sites. Call $`\kappa =k/N`$ the fraction of failed sites. Then from Eq. (3) the change in stress for any stable site $`i`$ is

$`\mathrm{\Delta }\sigma _{i\{k\}}`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{j\{k\}}{}}\mathrm{\Delta }_j={\displaystyle \frac{\delta }{N}}{\displaystyle \underset{j\{k\}}{}}(\sigma _j^f\sigma _j^R)`$ (5)
$`=`$ $`\delta \kappa \left(\sigma ^f_k\sigma ^R_k\right)`$ (6)
where $`\delta =(K_R+11/N)^1`$, and $`_k`$ is an average applied over failed sites. We call this the external stress transfer to signify that it applies to sites that are not part of the rupture. The quantity $`\sigma _j^f`$ is the stress of site $`j`$ at failure, which may be greater than $`\sigma ^F`$. Note that the slip displacement $`\mathrm{\Delta }_j`$ is only dependent on the stress drop at failure because pinning occurs immediately, and subsequent stress changes will not affect the slip. Since $`\sigma ^f`$ is typically very near $`\sigma ^F`$ and the $`\sigma _j^R`$ are identically distributed random variables, the term in parenthesis is on average independent of $`k`$. Thus the external transfer grows linearly (on average) with the fraction of failed sites.

In ‘dynamic weakening’ models  the increased propensity for a failed site to slip further (due to a weakened pinning force) is simulated by imposing a lower threshold stress $`\sigma ^D<\sigma ^F`$ for the duration of a single event. After a site fails it will receive stress transfer from subsequent failures, and thus may reach this lower threshold and fail again. Re-failing sites will contribute more to the external transfer and enhance the likelihood of continued rupture growth.

This form of weakening is strange in that it first requires some failed sites to have their stress brought back up to the dynamical threshold, which will occur at some minimum rupture size. Following the onset of dynamical weakening the additional stress transfer typically results in a runaway event which fails every block in the system. This feature can be exploited to produce characteristic events which always occur once the minimum size is reached.

One way to visualize the effects of weakening is to examine the average stress of sites that have failed as a function of rupture size (Figure 1). Without weakening, the average stress of failed sites (the average internal stress) will itself grow linearly like the external transfer (with half the slope). However, with dynamic weakening, all failed sites with stress $`\sigma ^D`$ will fail again, putting a ceiling on the average internal stress.

We claim that a more realistic approach would be to have failed sites shed a certain fraction of the stress they receive after failure. Implementing this ‘fractional weakening’ would involve re-computing the slip of all failed sites with each new failure. There is an easier way to accomplish this if we note that the desired effect is to lower the slope of the average internal stress (AIS) function. There is a way to invert this relationship, such that the average internal stress function is given and the requisite slips and stress transfer computed as a result. In essence, in place of solving dynamical equations involving slip or velocity dependent friction, we can impose the effects of the weakening as they appear in a discrete time context. We call this approach ‘forced weakening’.

To perform this inversion we first observe the change in stress of a site $`j`$ as it depends on $`k`$, the number of failed sites. Let $`\mathrm{\Delta }_j^k`$ denote the slip displacement of site $`j`$ after $`k`$ failures. If it is nonzero it includes all block motion, including initial failure, additional failures from weakening, or continuous sliding. Consider the system of equations for the stress changes of the failed sites $`j\{k\}`$

$`\mathrm{\Delta }\sigma _{j\{k\}}^k`$ $`=`$ $`\sigma _j^k\sigma _j^0=\delta \mathrm{\Delta }_j^k+{\displaystyle \frac{1}{N}}{\displaystyle \underset{ij\{k\}}{}}\mathrm{\Delta }_i^k`$ (7)
$`=`$ $`(K_R+1)\mathrm{\Delta }_j^k+{\displaystyle \frac{1}{N}}{\displaystyle \underset{i\{k\}}{}}\mathrm{\Delta }_i^k.`$ (8)
The last line demonstrates the simple linear form of the relationship between stress drops and slip displacements. This may be obtained via a matrix with diagonal elements $`N^1(K_R+1)`$ and off-diagonal elements $`N^1`$. This matrix may easily be inverted to obtain the slips $`\mathrm{\Delta }_j^k`$ in terms of the current stress drops $`\mathrm{\Delta }\sigma _j^k`$

$`\mathrm{\Delta }_j^k`$ $`=`$ $`{\displaystyle \frac{(K_R+1k/N)\mathrm{\Delta }\sigma _j^kN^1_{i\{k\}}\mathrm{\Delta }\sigma _i^k}{(K_R+1)\left(K_R+1k/N\right)}}`$ (9)
$`=`$ $`{\displaystyle \frac{\mathrm{\Delta }\sigma _j^k}{K_R+1}}{\displaystyle \frac{(k/N)\mathrm{\Delta }\sigma ^k_k}{(K_R+1)\left(K_R+1k/N\right)}}`$ (10)
This expression provides the slips $`\mathrm{\Delta }_j^k`$ that are necessary to generate a set of given stress drops.

The slips are directly related to the external transfer, as in Eq. 6. Summing over them yields a new expression for the external transfer in terms of the stress drops.

$`\mathrm{\Delta }\sigma _{i\{k\}}`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{j\{k\}}{}}\mathrm{\Delta }_j={\displaystyle \frac{k}{N}}\mathrm{\Delta }\sigma ^k_k`$ (12)
$`\times \left[{\displaystyle \frac{1}{K_R+1}}+{\displaystyle \frac{k/N}{(K_R+1)(K_R+1k/N)}}\right]`$
$`=`$ $`{\displaystyle \frac{\kappa }{(K_R+1\kappa )}}\mathrm{\Delta }\sigma ^k_\kappa `$ (13)
This is identical to Eq. (6) when $`k=1`$. However, the effective transfer coefficient $`\delta =\delta (\kappa )=(K_R+1\kappa )^1`$ now grows with the rupture size. Making up for this is the fact that the stress drops are no longer frozen after failure and will now decrease as the rupture grows. Observe that to calculate the external transfer, we need not specify individual stress drops or slips, but we just need to know the average dynamic stress drop of all failed sites for a given rupture size. This evolving average stress drop is defined as

$`\mathrm{\Delta }\sigma ^k_\kappa `$ $`=`$ $`{\displaystyle \frac{1}{k}}{\displaystyle \underset{j\{k\}}{}}(\sigma _j^0\sigma _j^k)`$ (14)
$`=`$ $`\sigma ^0_\kappa f(\kappa )`$ (15)
where we have defined the average internal stress (AIS) function $`f(\kappa )=\sigma ^k_k`$. The stress transfer now depends only on the initial stress and the AIS function, which are independent of failure order. Using this formulation, avalanches are Abelian.

In numerical simulation, we must eventually assign an actual slip and/or residual stress to each failed site. Care must be taken to make results agree with the analogous forward simulation. For example, given a linear AIS function $`f(\kappa )=\alpha \kappa +\beta `$, we could assign the $`k^{th}`$ failed site a random residual stress (with mean $`\beta `$) plus $`2\alpha (k1)/N`$. When using AIS functions with no forward equivalent, the method of assigning final stresses must be stated explicitly.

The forced weakening method has two main advantages. Practically, it allows the simulation of models with arbitrary weakening characteristics, some of which would not be obtainable with modified CA rules. Formally, the model is Abelian, so that the event size is a unique function of the initial stress configuration. Using this fact we can seek an expression which will determine the event size given adequate information of the initial stresses.

To accomplish this we first approximate the stress configuration with a continuous distribution $`p_\sigma (\chi )`$, where $`p_\sigma (\chi )d\chi `$ is the probability that a randomly selected site will have stress $`\sigma `$ between $`\chi `$, and $`\chi +d\chi `$. There are two ways of obtaining a continuous probability distribution from a discrete set of stresses: allow the number of sites $`N`$ to become infinite, or to consider averages over a course grained time. In the latter case the continuous driving of the slip deficits will quickly fill in the spectrum of possible stresses.

When considering a course grained time, we realize that the distribution is only an average, and that the actual instantaneous distribution of stresses will fluctuate. Similarly, the solution we obtain will be an average with statistical properties we would like to calculate. Under ordinary conditions, the distribution of stresses is statistically stationary , so the expected distribution of event sizes at any time will predict the long-time average behavior.

To write down the next-event-size expression, it is convenient to define the stress deficit $`\mathrm{\Sigma }=\sigma ^F\sigma `$, and the cumulative distribution

$$P_\mathrm{\Sigma }(\chi )=_0^\chi p_\mathrm{\Sigma }(\chi ^{})𝑑\chi ^{}$$
(16)
where $`p_\mathrm{\Sigma }(\chi )=p_\sigma (\sigma ^F\chi )`$. Given this distribution the event size $`\kappa `$ is the solution to

$`{\displaystyle \frac{1}{K_R+1\kappa }}\left[\kappa \sigma ^F\kappa f(\kappa )\sigma ^F{\displaystyle _0^\kappa }P_\mathrm{\Sigma }(\chi )𝑑\chi \right]`$ (17)
$`\sigma ^FP_\mathrm{\Sigma }^1(\kappa +d\kappa )`$ $`=`$ $`0`$ (18)
where $`f(\kappa )`$ is again the AIS function. A full derivation and analysis of the results of this equation will be published elsewhere. For an example solution, suppose each stress in the system is chosen independently from a uniform distribution between zero and $`\sigma ^F`$. Then for large $`N`$, Eq. (18) determines that $`\kappa `$ is the first crossing of a 1-d random walk constrained to return to zero at $`\kappa =1`$. Without that constraint, this produces a power law distribution of event sizes with exponent $`3/2`$. With the constraint in place, there is an upturn at large event sizes, which is observed in simulations (so long as $`K_R`$ is low enough to allow such large events). A comparison of this analytical solution with a simulation is shown in Figure 2.

In summary, the forced weakening method extends the basic slider block model to include arbitrary weakening functions in an efficient discrete time simulation. Additionally, forced weakening results in Abelian avalanches which leads to a rigorous and implementation-independent analysis. We have presented a brief example of this showing excellent agreement between theory and simulation. This approach should be equally applicable to related discrete threshold models. While the resulting formalism is dependent on the mean-field character of the model, the numerical techniques may find wider use whenever an inversion like Eq. (10) is available.

The authors thank W. Klein for his input. E.F.P. and J.B.R. acknowledge support by DOE grant DE-FG03-95ER14499; J.S.S.M. was supported as a Visiting Fellow by CIRES, University of Colorado at Boulder.


text:
# Testing Lorentz and CPT symmetry with hydrogen masers

## I Introduction

A theoretical framework has recently been developed that incorporates Lorentz and CPT symmetry violation into the standard model and quantifies their effects. . One branch of this framework emphasizes low energy, experimental searches for symmetry violating effects in atomic energy levels . In particular, Lorentz and CPT violation in hydrogen has been examined and sidereal variations in the $`F=1`$, $`\mathrm{\Delta }m_F=\pm 1`$ Zeeman frequency have been quantified . Motivated by this work, we have conducted a search for sidereal variation in the hydrogen Zeeman frequency, and have placed a new clean bound of 10-27 GeV on Lorentz and CPT violation of the proton .

Here we provide additional details of the theoretical framework, experiment and analysis. In Sec. II we discuss the standard model extension. In Sec. III we describe the basic concepts of hydrogen maser operation and our Zeeman frequency measurement technique. In Sec. IV we describe the procedure used to collect data and extract a sidereal bound on the Zeeman frequency. In Sec. V we describe efforts to reduce and characterize systematic effects. Finally, in Sec. VI we compare our result to other clock-comparison tests of Lorentz and CPT symmetry, and discuss potential means of improving our measurement.

## II Lorentz and CPT symmetry violation in the standard model

Experimental investigations of Lorentz symmetry provide important tests of the standard model of particle physics and general relativity. While the standard model successfully describes particle phenomenology, it is believed to be the low energy limit of a fundamental theory that incorporates gravity. This underlying theory may be Lorentz invariant, yet contain spontaneous symmetry-breaking that could result in small violations of Lorentz invariance and CPT at the level of the standard model.

A theoretical framework has been developed to describe Lorentz and CPT violation at the level of the standard model by Kostelecký and coworkers . This standard-model extension is quite general: it emerges as the low-energy limit of any underlying theory that generates the standard model and contains spontaneous Lorentz symmetry violation . For example, such characteristics might emerge from string theory . A key feature of the standard model extension is that it is formulated at the level of the known elementary particles, and thus enables quantitative comparison of a wide array of searches for Lorentz and CPT violation .

“Clock comparison experiments” are searches for temporal variations in atomic energy levels. According to the standard model extension considered here, Lorentz and CPT violation may produce shifts in certain atomic levels, whose magnitude depends on the orientation of the atom’s quantization axis relative to a fixed inertial frame . Certain atomic transition frequencies, therefore, may exhibit sinusoidal variation as the earth rotates on its axis. New limits can be placed on Lorentz and CPT violation by bounding sidereal variation of these atomic transition frequencies.

Specifically, the description of Lorentz and CPT violation is included in the relativistic Lagrange density of the constituent particles of the atom. For example, the modified electron Lagrangian becomes

$$=\frac{1}{2}i\overline{\psi }\mathrm{\Gamma }_\nu ^\nu \psi \overline{\psi }M\psi +_{int}^{QED}$$
(1)
where

$$\mathrm{\Gamma }_\nu =\gamma _\nu +\left(c_{\mu \nu }\gamma ^\mu +d_{\mu \nu }\gamma _5\gamma ^\mu +e_\nu +if_\nu \gamma _5+\frac{1}{2}g_{\lambda \mu \nu }\sigma ^{\lambda \mu }\right)$$
(2)
and

$$M=m+\left(a_\mu \gamma ^\mu +b_\mu \gamma _5\gamma ^\mu +\frac{1}{2}H_{\mu \nu }\sigma ^{\mu \nu }\right).$$
(3)
The parameters $`a_\mu `$, $`b_\mu `$, $`c_{\mu \nu }`$, $`d_{\mu \nu }`$, $`e_\nu `$, $`f_\nu `$, $`g_{\lambda \mu \nu }`$ and $`H_{\mu \nu }`$ represent possible vacuum expectation values of Lorentz tensors generated through spontaneous Lorentz symmetry breaking in an underlying theory. These are absent in the standard model. The parameters $`a_\mu `$, $`b_\mu `$, $`e_\nu `$, $`f_\nu `$ and $`g_{\lambda \mu \nu }`$ represent coupling strengths for terms that violate both CPT and Lorentz symmetry, while $`c_{\mu \nu }`$, $`d_{\mu \nu }`$, and $`H_{\mu \nu }`$ violate Lorentz symmetry only. An analogous expression exists for the modified proton and neutron Lagrangians (a superscript will be appended to differentiate between the sets of parameters). The standard model extension treats only the free particle properties of the constituent particles, estimating that all interaction effects will be of higher order . As a result, the interaction term $`_{int}^{QED}`$ is unchanged from the conventional, Lorentz invariant, QED interaction term.

Within this phenomenological framework, the values of these parameters are not calculable; instead, values must be determined experimentally. The general nature of this theory ensures that different experimental searches may place bounds on different combinations of Lorentz and CPT violating terms, while direct comparisons between these experiments are possible (see Table III and Ref. ).

The leading-order Lorentz and CPT violating energy level shifts for a given atom are obtained by summing over the individual free particle shifts of the atomic constituents. From the symmetry violating correction to the relativistic Lagrangian, a non-relativistic correction Hamiltonian $`\delta h`$ is found using standard field theory techniques . Assuming Lorentz and CPT violating effects to be small, the energy level shifts are calculated perturbatively by taking the expectation value of the correction Hamiltonian with respect to the unperturbed atomic states, leading to a shift in an atomic ($`F,m_F`$) sublevel given by

$$\mathrm{\Delta }E_{F,m_F}=F,m_F|n_e\delta h_e+n_p\delta h_p+n_n\delta h_n|F,m_F.$$
(4)
Here $`n_w`$ is the number of each type of particle and $`\delta h_w`$ is the corresponding correction Hamiltonian. Note that for most atoms, the interpretation of energy level shifts in terms of this standard model extension is reliant on the particular model used to describe the atomic nucleus (e.g., the Schmidt model). One key advantage of a study in hydrogen is the simplicity of the nuclear structure (a single proton), with its results uncompromised by any nuclear model uncertainty.

Among the most recent clock comparison experiments are Penning trap tests by Dehmelt and co-workers with the electron and positron  which place a limit on electron Lorentz and CPT violation at 10-25 GeV. A recent re-analysis by Adelberger, Gundlach, Heckel, and co-workers of existing data from the “Eöt-Wash II” spin-polarized torsion pendulum  has improved this to a level of 10-29 GeV , the most stringent bound to date on Lorentz and CPT violation of the electron. A new limit on neutron Lorentz and CPT violation has been placed at 10-31 GeV by Bear et al.  using a dual species noble gas maser and comparing Zeeman frequencies of 129Xe and 3He. The current limit on Lorentz and CPT violation of the proton is 10-27 GeV, as derived from an experiment by Lamoreaux and Hunter  which compared Zeeman frequencies of 199Hg and 133Cs.

Figure 1 shows the Lorentz and CPT violating corrections to the energy levels of the ground state of hydrogen . The shift in the $`F=1`$, $`\mathrm{\Delta }m_F=\pm 1`$ Zeeman frequency is :

$$|\mathrm{\Delta }\nu _Z|=\frac{1}{h}|(b_3^ed_{30}^em_eH_{12}^e)+(b_3^pd_{30}^pm_pH_{12}^p)|.$$
(5)
The subscripts denote the projection of the tensor couplings onto the laboratory frame. Therefore, as the earth rotates relative to a fixed inertial frame, the Zeeman frequency $`\nu _Z`$ will exhibit a sidereal variation. We have recently published the result of a search for this variation of the $`F=1`$, $`\mathrm{\Delta }m_F=\pm 1`$ Zeeman frequency in hydrogen using hydrogen masers . This search has placed a new, clean bound on Lorentz and CPT violation of the proton at a level of 10-27 GeV.

## III Hydrogen maser concepts

The electronic ground state in hydrogen is split into four levels by the hyperfine interaction, labeled (following the notation of Andresen ) $`|1`$ to $`|4`$ in order of decreasing energy (Fig. 1). The energies of atoms in $`|1`$ and $`|2`$ decrease as the magnetic field decreases; these are therefore low-field seeking states. Conversely, $`|3`$ and $`|4`$ are high-field seeking states. In low fields, $`|2`$ and $`|4`$ are only dependent on magnetic field in second order. The maser oscillates on the $`|2|4`$ transition (field-independent to first-order). This transition frequency, as a function of static field, is given by $`\nu _{24}=\nu _{hfs}+2750B^2`$ ($`\nu `$ in Hz with $`B`$ in Gauss, with $`\nu _{hfs}`$ 1420.405751 MHz the zero-field hyperfine frequency). Hydrogen masers typically operate with low static fields (less than 1 mG), where $`\nu _{24}`$ is shifted from $`\nu _{hfs}`$ by about 3 mHz, or 2 parts in 1012. The two $`F=1`$, $`\mathrm{\Delta }m_F=\pm 1`$ Zeeman frequencies are given by $`\nu _{12}=1.4\times 10^6B1375B^2`$ and $`\nu _{23}=1.4\times 10^6B+1375B^2`$. At $`B`$ = 1 mG these are nearly degenerate, with $`\nu _{12}\nu _{23}`$ 3 mHz, much less than the Zeeman linewidth of approximately 1 Hz.

### A Maser operation

In a hydrogen maser , molecular hydrogen is dissociated in an rf discharge and a beam of hydrogen atoms is formed, as shown in Fig. 2. A hexapole state selecting magnet focuses the low-field-seeking hyperfine states $`|1`$ and $`|2`$ into a quartz maser bulb at about $`10^{12}`$ atoms/sec. Inside the bulb (volume $`10^3`$ cm3), the atoms travel ballistically for about 1 second before escaping, making $`10^4`$ collisions with the bulb wall. A Teflon coating reduces the atom-wall interaction and thus inhibits decoherence of the masing atomic ensemble by wall collisions. The maser bulb is centered inside a cylindrical TE011 microwave cavity resonant with the 1420 MHz hyperfine transition. The microwave field stimulates a small, coherent magnetization in the atomic ensemble, and this magnetization acts as a source to stimulate the microwave field. With sufficiently high atomic flux and low cavity losses, this feedback induces active maser oscillation. The maser signal is inductively coupled out of the microwave cavity and amplified with an external receiver. Surrounding the cavity, a solenoid produces the weak static magnetic field ($``$ 1 mG) that establishes the quantization axis inside the maser bulb and sets the Zeeman frequency ($``$ 1 kHz). A pair of Helmholtz coils produces the oscillating transverse magnetic field that drives the $`F=1`$, $`\mathrm{\Delta }m_F=\pm 1`$ Zeeman transitions. The cavity, solenoid, and Zeeman coils are all enclosed within several layers of high permeability magnetic shielding.

A well engineered hydrogen maser can have fractional stabilities approaching $`10^{15}`$ over intervals of hours. This stability is enabled by a long atom-field interaction time (1 s), a low atom-wall interaction (due to the low atomic polarizability of H and the wall’s Teflon coating), reduced Doppler effects (the atoms are confined to a region of uniform microwave field phase, effectively averaging their velocity to zero over the interaction time with the field), and multiple layers of thermal control of the cavity (stabilizing cavity pulling shifts).

### B Maser characterization

Among the quantities used to characterize a hydrogen maser, those most relevant to this experiment are the atomic line-Q $`Q_l`$, the population decay rate $`\gamma _1`$, the hyperfine decoherence rate $`\gamma _2`$, the atomic flow rate into and out of the bulb $`\gamma _b`$, and the maser Rabi frequency $`|X_{24}|`$. We describe here a comprehensive set of measurements to characterize hydrogen maser P-8. The results discussed here are summarized in Table I. Our Lorentz and CPT symmetry test data were taken with a similar but newer hydrogen maser, P-28 . A few of the maser characterization parameters for P-28, while not directly measured, have been inferred using fitting parameters from the double resonance method used to measure the $`F=1`$, $`\mathrm{\Delta }m_F=\pm 1`$ Zeeman frequency, described in Sec. IV A. These values are included in Table I in italics.

To determine these parameters of an operating H maser, the cavity volume $`V_C`$, bulb volume $`V_b`$, cavity quality factor $`Q_C`$, filling factor $`\eta `$, and output coupling coefficient $`\beta `$ must be known. For both masers, $`V_C`$ = 1.4 $`\times `$ 10-2 m3, $`V_b`$ = 2.9 $`\times `$ 10-3 m3, $`Q_C`$ 40,000, and $`\beta `$ = 0.23 . The filling factor, defined as

$$\eta =\frac{H_z_{bulb}^2}{H^2_{cavity}},$$
(6)
quantifies the ratio of average magnetic field energy inside the bulb to the average magnetic field energy in the cavity. This has a value of $`\eta `$ = 2.14 for masers P-8 and P-28 .

For canonical hydrogen maser operation, there are two important relaxation rates . For a room temperature H maser, the decay of the population inversion is described by the longitudinal relaxation rate

$$\gamma _1=\gamma _b+\gamma _r+2\gamma _{se}+\gamma _1^{},$$
(7)
and the decay of the atomic coherence is described by the transverse relaxation rate

$$\gamma _2=\gamma _b+\gamma _r+\gamma _{se}+\gamma _2^{}.$$
(8)
Here, $`\gamma _b`$ is the atomic flow rate into the bulb, $`\gamma _r`$ is the rate of recombination into molecular hydrogen at the bulb wall, $`\gamma _{se}`$ is the hydrogen-hydrogen spin-exchange decay rate, and $`\gamma _i^{}`$ includes all other sources of decay, such as decoherence during wall collisions and effects of magnetic field gradients.

In the steady state, the atom flow rate into the bulb is equal to the geometric escape rate from the bulb, given by $`\gamma _b=\overline{v}A/4KV_b`$, where $`\overline{v}`$ = 2.5 $`\times 10^5`$ cm/s is the mean thermal velocity of atoms in the bulb, $`A`$ = 0.254 cm3 is the area of the bulb entrance aperture, and $`K`$ 6 is the Klausing factor . Thus, $`\gamma _b`$ = 0.86 rad/s for both P-8 and P-28. The spin exchange decay rate is given approximately by

$$\gamma _{se}=\frac{1}{2}n\overline{v}_r\sigma $$
(9)
where $`\overline{v}_r`$ = 3.6 $`\times 10^5`$ cm/s is the mean relative velocity of atoms in the bulb and $`\sigma `$ = 21 $`\times `$ 10-16 cm2 is the hydrogen-hydrogen spin-exchange cross section. The hydrogen density is given by

$$n=\frac{I_{tot}}{(\gamma _b+\gamma _r)V_b}$$
(10)
where $`I_{tot}`$ is the total flux of hydrogen atoms into the storage bulb.

The atomic line-Q is related to the transverse relaxation rate and the maser oscillation frequency $`\omega `$ by

$$Q_l=\frac{\omega }{2\gamma _2}.$$
(11)
It is measured using the cavity pulling of the maser frequency: neglecting spin-exchange shifts, the maser frequency is given by

$$\omega =\omega _{24}+\frac{Q_C}{Q_l}\left(\omega _C\omega _{24}\right).$$
(12)
By measuring the maser frequency as a function of cavity frequency setting, the line-Q can be determined. For both P-8 and P-28, we find $`Q_l`$ = 1.6 $`\times `$ 109, and therefore $`\gamma _2`$ = 2.8 rad/s.

A convenient single measure of spin-exchange-independent relaxation in a hydrogen maser is given by “gamma-t”

$$\gamma _t=\left[(\gamma _b+\gamma _r+\gamma _1^{})(\gamma _b+\gamma _r+\gamma _2^{})\right]^{\frac{1}{2}}.$$
(13)
Using this, a more useful form for the longitudinal relaxation rate, $`\gamma _1`$, can be found. By combining Eqn. 13 with Eqns. 7 and 8, we find

$$\gamma _1=\frac{\gamma _t^2}{\gamma _2\gamma _{se}}+2\gamma _{se}.$$
(14)
Using Eqns. 8-11, we can relate the line-Q to $`I`$, the input flux of atoms in state $`|2`$ as

$$\frac{1}{Q_l}=\frac{2}{\omega }\left[\gamma _b+\gamma _r+\gamma _2^{}+q\frac{I}{I_{th}}\gamma _t\right]$$
(15)
using the threshold flux required for maser oscillation (neglecting spin-exchange)

$$I_{th}=\frac{\mathrm{}V_C\gamma _t^2}{4\pi \mu _B^2Q_C\eta },$$
(16)
and the maser quality parameter

$$q=\left[\frac{\sigma \overline{v}_r\mathrm{}}{8\pi \mu _b^2}\right]\frac{\gamma _t}{\gamma _b+\gamma _r}\left[\frac{V_C}{\eta V_b}\right]\left(\frac{1}{Q_C}\right)\frac{I_{tot}}{I}.$$
(17)
The ratio $`I`$/$`I_{tot}`$ is a measure of the effectiveness of the state selection of atoms entering the bulb. While $`I`$ is not directly measurable, it can be related to the power $`P`$ radiated by the atoms by

$$\frac{P}{P_c}=2q^2\left(\frac{I}{I_{th}}\right)^2+(13q)\left(\frac{I}{I_{th}}\right)1$$
(18)
where $`P_c=\mathrm{}\omega I_{th}/2`$. The maser power is also related to the maser Rabi frequency by

$$P=\frac{I\mathrm{}\omega }{2}\frac{|X_{24}|^2}{\gamma _1\gamma _2}\left(1+\frac{|X_{24}|^2}{\gamma _1\gamma _2}\right)^1.$$
(19)
The power coupled out of the maser is given by  $`P_o/P=\beta /(1+\beta )`$.

Generally, the parameter $`q`$ is less than 0.1, while $`I/I_{th}`$ is approximately 2 or 3. Hence, the first term of Eqn. 18 can be neglected relative to the others. If we make the reasonable approximation that $`\gamma _1^{}`$ = $`\gamma _2^{}`$, then we can rewrite Eqn. 15 using Eqns. 13, 16-18 as

$$\frac{1}{Q_l}=mP+b$$
(20)
where

$$b=\frac{2}{\omega }\gamma _t\left[1+\frac{q}{13q}\right]$$
(21)
and

$$m=\frac{16\pi \mu _b^2Q_C\eta }{\omega ^2\mathrm{}^2V_C}\left[\frac{q}{13q}\right]\frac{1}{\gamma _t}.$$
(22)
Therefore, by measuring the line-Q as a function of maser power and extracting the slope $`m`$ and the y-intercept $`b`$, we can determine $`q`$ and $`\gamma _t`$. For maser P-8, $`q`$ = 0.100 and $`\gamma _t`$ = 1.83 rad/s.

With these values of $`q`$ and $`\gamma _t`$, we found $`I_{th}`$ = 0.54 $`\times `$ 1012 atoms/s (using Eqn. 16), and $`P_c`$ = 250 fW. With a measured output power of $`P_o`$ = 112 fW, the atoms were radiating $`P`$ = 599 fW, and the flux of state $`|2`$ atoms was $`I`$ = 3.13 $`\times `$ 1012 atoms/s (Eqn. 18). Under the assumption that $`\gamma _t\gamma _b+\gamma _r`$ we found that the total flux was $`I_{tot}`$ = 15.0 $`\times `$ 1012 atoms/s (Eqn. 17) and the density was $`n`$ = 2.8 $`\times `$ 1015 atoms/m3 (Eqn. 10). The spin-exchange decay rate was then found to be $`\gamma _{se}`$ = 1.06 rad/s (Eqn. 9). Finally, the population decay rate was $`\gamma _1`$ = 4.04 rad/s (Eqn. 14) and the maser Rabi frequency was $`|X_{24}|`$ = 2.77 rad/s (Eqn. 19).

### C Zeeman frequency determination

The $`F=1`$, $`\mathrm{\Delta }m_F=\pm 1`$ Zeeman frequency is measured using a double resonance technique . As the frequency of an audio frequency magnetic field $`\omega _Z`$, applied perpendicular to the quantization axis, is swept through the Zeeman frequency, a shift in the maser frequency is observed (Fig. 3). When the applied field is near the Zeeman frequency, two-photon transitions (one audio photon plus one microwave photon) link states $`|1`$ and $`|3`$ to state $`|4`$, in addition to the single microwave photon transition between states $`|2`$ and $`|4`$. This two photon coupling shifts the maser frequency antisymmetrically with respect to the detuning of the applied field about the Zeeman resonance .

To second order in the Rabi frequency of the applied Zeeman field, $`|X_{12}|`$ the small static-field limit of the maser frequency shift from the unperturbed frequency is given by

$`\mathrm{\Delta }\omega `$ $`=`$ $`|X_{12}|^2(\rho _{11}^0\rho _{33}^0){\displaystyle \frac{\delta (\gamma _1\gamma _2+|X_{24}^0|^2)(\gamma _Z/\gamma _b)}{(\gamma _Z^2\delta ^2+\frac{1}{4}|X_{24}^0|^2)^2+(2\delta \gamma _Z)^2}}`$ (24)
$`+|X_{12}|^2\left({\displaystyle \frac{\omega _C\omega _{24}}{\omega _{24}}}\right){\displaystyle \frac{Q_C\gamma _Z(1+K)}{\gamma _Z^2(1+K)^2+\delta ^2(1K)^2}}`$
where $`\gamma _Z`$ is the Zeeman decoherence rate, $`\delta =\omega _Z\omega _{23}`$ is the detuning of the applied field from the atomic Zeeman frequency, $`K=\frac{1}{4}|X_{24}^0|^2/(\gamma _Z^2+\delta ^2)`$, and $`\rho _{11}^0\rho _{33}^0=\gamma _b/(2\gamma _1)`$ is the steady state population difference between states $`|1`$ and $`|3`$ in the absence of the applied Zeeman field. The first term in Eqn. 24 results from the coherent two-photon mixing of the $`F=1`$ levels as described above , while the second term is a modified cavity pulling term that results from the reduced line-Q in the presence of the applied Zeeman field. We compared Eqn. 24 to experimental data from P-8, inserting the independently measured values of $`|X_{24}^0|`$, $`\gamma _b`$, $`\gamma _1`$, and $`\gamma _2`$. By matching the fit to the data we extracted the Zeeman field parameters $`|X_{12}|`$ and $`\gamma _Z`$ shown in Fig. 3.

In addition to the shift given by Eqn. 24, there is a small symmetric frequency shift due to the slight non-degeneracy of the two $`F=1`$, $`\mathrm{\Delta }m_F=\pm 1`$ Zeeman frequencies. This term offsets the zero crossing of the maser shift resonance away from the average Zeeman frequency $`\frac{1}{2}\left(\nu _{12}+\nu _{23}\right)`$, however the contribution is negligible at small static fields. Also, a reanalysis of the double resonance maser shift , which included the effects of spin-exchange collisions , showed that there is an additional hydrogen density-dependent offset of the zero crossing of the maser shift resonance from the average Zeeman frequency. Using the full spin-exchange corrected formula for the maser frequency shift , we calculated this offset and found that for typical hydrogen maser densities ($`n`$ 3 $`\times 10^{15}`$ m-3), the offset varied with average maser power as approximately -50 $`\mu `$Hz/fW (assuming a linear relation between maser power and atomic density of $`\frac{\mathrm{\Delta }P}{\mathrm{\Delta }n}\frac{100fW}{3\times 10^{15}m^3}`$). As described below, our masers typically have sidereal power fluctuations less than 1 fW, making this effect negligible.

The applied Zeeman field also acts to diminish the maser power, as shown in Fig. 4, and to decrease the maser’s line-Q. By driving the $`F=1`$, $`\mathrm{\Delta }m_F=\pm 1`$ Zeeman transitions, the applied field depletes the population of the upper masing state $`|2`$, thereby diminishing the number of atoms undergoing the maser transition and reducing the maser power. Also, by decreasing the lifetime of atoms in state $`|2`$, the line-Q is reduced. A very weak Zeeman field of about 35 nG (as was used in our Lorentz symmetry test) decreases the maser power by less than 2% on resonance and reduces the line-Q by 2% (as calculated using Eqn. 6 of ). The standard method of determining the average static magnetic field strength is to scan the Zeeman resonance with a large applied field and record the power diminishment (such as that shown in Fig. 4, open circles). From the applied field frequency at the center of the power resonance, which typically has a width of about 1 Hz, the magnetic field can be found with a resolution of about 1 $`\mu `$G.

## IV Experimental procedure

### A Zeeman frequency measurement

To measure the $`F=1`$, $`\mathrm{\Delta }m_F=\pm 1`$ Zeeman frequency, we applied an oscillating field of about 35 nG near the Zeeman frequency. This field shifted the maser frequency by a few mHz (at the extrema), a fractional shift of about 2 parts per trillion. Because of the excellent fractional maser stability (2 parts in $`10^{14}`$ over our averaging times of 10 s), the shift was easily resolved (see the solid data in Fig. 3). As the frequency of the applied field was stepped through the Zeeman resonance, the maser frequency (of perturbed maser P-28) was compared to a second, unperturbed hydrogen maser frequency (P-13). The two maser signals at $``$ 1420 MHz were phase locked to independent voltage controlled crystal oscillator receivers. The exact value of the receivers’ outputs were set by tunable synthesizers, which were set such that there was a 1.2 Hz offset between them. The two receiver outputs were combined in a heterodyne mixer and the resulting 0.8 s period beat note was averaged for 10 s (about 12 periods) with a Hewlett-Packard Model HP 5334B frequency counter. The full double resonance spectrum consisted of 100 such points. For each spectrum, 80% of the points were taken over the middle 40% of the scan range, where the frequency shift varies the most.

Once an entire spectrum of beat period vs applied Zeeman frequency was obtained, it was fit to the function

$`T_b`$ $`=A_0`$ $`+{\displaystyle \frac{A_3\delta (1\kappa )}{A_1(1+\kappa )^2+\delta ^2(1\kappa )^2}}{\displaystyle \frac{A_3(\delta +\tau )(1\kappa )}{A_1(1+\kappa )^2+(\delta +\tau )^2(1\kappa )^2}}`$ (26)
$`+{\displaystyle \frac{A_5\delta }{(A_1\delta ^2+A_4)^2+4\delta ^2A_1}}+{\displaystyle \frac{A_6(1+\kappa )}{A_1(1+\kappa )^2+\delta ^2(1\kappa )^2}}`$
to determine the Zeeman frequency. Here $`\delta =\nu \nu _Z`$ is the Zeeman detuning of the applied field $`\nu `$ away from the Zeeman frequency $`\nu _Z`$, $`\kappa =A_4/(A_1+\delta ^2)`$ is the analog of the parameter $`K`$ from Eqn. 24, and $`\tau =(1.403\times 10^9)\times \nu _Z^2`$ is the small difference between the two Zeeman frequencies $`\nu _{12}`$ and $`\nu _{23}`$. The first term $`A_0`$ is the constant offset representing the unperturbed beat period between the two masers. The second and third terms comprise the first-order symmetric maser shift (not included in Eqn. 24 but described in the text above); these two terms nearly cancel at low static field where $`\tau `$ vanishes. The final two terms account for the two shifts given in Eqn. 24.

For our spectra in maser P-28 with small applied field amplitude (solid square data points of Fig. 3), typical fit parameters were: $`A_0`$ = 0.84550 $`\pm `$ 0.00001, $`A_1`$ = 0.141 $`\pm `$ 0.005, $`\nu _Z`$ = 857.063 $`\pm `$ 0.003, $`A_3`$ = 0.006 $`\pm `$ 0.010, $`A_4`$ = 0.029 $`\pm `$ 0.003, $`A_5`$ = (3.2 $`\pm `$ 0.1) $`\times 10^4`$, and $`A_6`$ = (-1 $`\pm `$ 5) $`\times 10^6`$. The uncertainty in the Zeeman frequency was 3 mHz. Also, $`A_3`$ and $`A_6`$, the amplitude coefficients of the residual first order effect and the cavity pulling term, were consistent with zero.

With our known value of $`\gamma _b`$ = 0.86 rad/s, and our measured value of $`\gamma _2`$ = 2.77 rad/s (from line-Q), the above set of fit parameters were consistent with the reasonable values $`X_{24}^0`$ = 2.14 rad/s, $`\gamma _Z`$ = 2.36 rad/s, $`\gamma _1`$ = 2.88 rad/s, and $`X_{12}`$ = 0.40 rad/s (since $`A_3`$ had such a large error bar, the value of $`X_{12}`$ was chosen such the ratio of the maser shift amplitude in P-28 to P-8, shown in Fig. 3, is equal to the ratio of the squares of $`X_{12}`$ for P-28 to P-8).

To determine the number of points and length of averaging that optimized the Zeeman frequency resolution, we recorded several spectra with 50, 100, and 150 points at 5 s and 10 s averaging. We also varied the “density distribution” of points, including spectra where the middle 40% of the scan contained 80% of the points and those where the middle 30% contained 80% of the points (thus increasing the number of points in the region where the antisymmetric shift varies the most). With each of these spectra, we ran the following Monte Carlo analysis : after fitting each scan to Eqn. 26, we constructed 100 synthetic data sets by adding Gaussian noise to the fit, with noise amplitude determined by the unperturbed maser frequency resolution of about 40 $`\mu `$Hz. Each of these synthetic data sets was fit and a histogram of the fitted Zeeman frequencies was constructed. The resolution of each spectrum was taken as the width of the Gaussian curve that fit the histogram (see Fig. 5). As the total length of the scans increased, the resolution improved and converged to a limit of around 2.5 mHz. While the resolution improved slowly with increased acquisition time, it would have eventually begun to degrade due to long term drifting of the Zeeman frequency. (As will be described below, we found that the Zeeman frequency exhibited slow drifts of about 10-100 mHz/day). We therefore chose a scan of 100 points at 10 s averaging, for a total length of about 18 minutes for our Lorentz symmetry test spectra. The results from the Monte Carlo analysis for one of these spectra indicated a Zeeman frequency resolution is 2.7 mHz (see Fig. 5).

### B Data analysis

Our net result combines data from three runs. During each data run, the 18 minute Zeeman frequency scans were automated and run consecutively. After every 10 scans, 20 minutes of “unperturbed” maser frequency stability data was taken to track the maser’s stability. Each run contained about 10 continuous days worth of data, and each set contained more than 500 Zeeman frequency measurements, taken at $``$ 18 minute intervals.

For each run, the long term Zeeman frequency data was fit to a function of the form

$$fit=(\mathrm{𝑝𝑖𝑒𝑐𝑒𝑤𝑖𝑠𝑒}\mathrm{𝑐𝑜𝑛𝑡𝑖𝑛𝑢𝑜𝑢𝑠}\mathrm{𝑙𝑖𝑛𝑒𝑎𝑟}\mathrm{𝑓𝑢𝑛𝑐𝑡𝑖𝑜𝑛})+\delta \nu _{Z,\alpha }\mathrm{cos}(\omega _{sid}t)+\delta \nu _{Z,\beta }\mathrm{sin}(\omega _{sid}t)$$
(27)
where $`\delta \nu _{Z,\alpha }`$ and $`\delta \nu _{Z,\beta }`$ represent the cosine and sine components of the sidereal sinusoid. The time origin of the sinusoids for all three runs was taken as midnight (00:00) of November 19, 1999. The subscripts $`\alpha `$ and $`\beta `$ refer to two non-rotating orthogonal axes perpendicular to the rotation axis of the earth. The total sidereal amplitude was determined by adding $`\delta \nu _{Z,\alpha }`$ and $`\delta \nu _{Z,\beta }`$ in quadrature. During each run, the Zeeman frequency drifted hundreds of mHz over tens of days. The piecewise continuous linear function, consisting of segments one sidereal day in length, was included to account for these long term Zeeman frequency drifts. This function was continuous at each break, while the derivative was discontinuous.

The result of this analysis, where the fitting function (Eqn. 27) was applied to the full data set, was found to be in good agreement with a second analysis, where each individual day of data was fit to a line plus the sidereal sinusoid and the cosine and sine amplitudes of each day were averaged separately and then combined in quadrature to find the total sidereal amplitude.

### C Run 1

The cumulative data from the first run (November 1999) are shown in Fig. 6(a) and the residuals from the complete fit (Eqn. 27) are shown in Fig. 6(b). The data set consisted of 11 full days of data and had an overall drift of about 250 mHz.

To avoid a biased choice of fitting, we allowed the location of the slope discontinuities in the piecewise continuous linear function to shift throughout a sidereal day. We made eight separate fits, each with the location of the slope discontinuities shifted by three sidereal hours. The total sidereal amplitude and reduced chi square for each is shown in Fig. 7. We chose our result from the fit with minimum reduced chi square.

As noted above, the error bar on a single Zeeman frequency determination was about 3 mHz. However, when analyzing a smooth region of long term Zeeman data (about 1 day) we calculate a standard deviation of about 5 mHz. We believe this error bar is due mainly to residual thermal fluctuations (see Fig. 14).

For our choice of slope discontinuity with minimum reduced chi square , the cosine amplitude was 0.43 mHz $`\pm `$ 0.36 mHz, and the sine amplitude was -0.21 mHz $`\pm `$ 0.36 mHz. The total sidereal amplitude was therefore 0.48 mHz $`\pm `$ 0.36 mHz.

### D Field-inverted runs 2 and 3

In runs 2 and 3, the static solenoid field orientation was opposite that of the initial run to further study the double resonance technique and any potential systematics associated with the solenoid field. With the static field inverted, and therefore directed opposite the quantization axis in the state selecting hexapole magnet, the input flux consists of atoms in states $`|2`$ and $`|3`$ (rather than the states $`|1`$ and $`|2`$). Thus, reversing the field inverts the steady state population difference $`(\rho _{11}^0\rho _{33}^0)`$ of Eqn. 24 and acts to invert the antisymmetric double resonance maser frequency shift .

Operating the maser in the field reversed mode degrades the maser performance and subsequently the Zeeman frequency data. With opposed quantization fields inside the maser bulb and at the exit of the state selecting hexapole magnet, a narrow region of field inversion is created. Where the field passes through zero, Majorana transitions between the different $`m_F`$ sublevels of the $`F=1`$ manifold can occur. This can alter the number of atoms in the upper maser state ($`F=1`$, $`m_F=0`$, state $`|2`$), which diminishes the overall maser amplitude and stability. In the field-inverted configuration, the maser amplitude was reduced by 30%, and both the maser frequency and Zeeman frequency were less stable. In addition, the field-inverted runs were each conducted soon after a number of rather invasive repairs were made to the maser . Thus, the quality of the latter two data sets was somewhat degraded from the first run (see Figs. 8(a) and 9(a)). The overall drift was larger (nearly 800 mHz over about 10 days), and the scatter in the data was increased, as can be seen from the residual plots from these runs (Figs. 8(b) and 9(b)) which have been plotted on the same scale as the residuals from the first run (Fig. 6(b)).

The latter two runs were also less suitable for the piecewise continuous linear drift model used in the first run. In that case, the large slope changes were coincidentally separated by an integer number of sidereal days; in the last two runs, the larger and more frequent changes in slope were not. Therefore, only certain selected sections could be fit to the same model (Eqn. 27), significantly truncating the data sets. Due to all of these factors, the sidereal amplitudes and the associated error bars were up to an order of magnitude larger for the field-inverted runs than the first run. All values are shown together in Table II.

### E Combined result

The final sidereal bound, combining all three runs, was calculated using the data in Table II. First, the weighted averages of the cosine and sine amplitudes, $`\overline{\delta \nu }_{Z,\alpha }`$ and $`\overline{\delta \nu }_{Z,\beta }`$, were found using the standard formula for weighted mean

$$\mu ^{}=\left(\mathrm{\Sigma }\frac{x_i}{\sigma _i^2}\right)/\left(\mathrm{\Sigma }\frac{1}{\sigma _i^2}\right),$$
(28)
and their uncertainties were given by

$$\sigma _\mu ^{}^2=\sqrt{1/\left(\mathrm{\Sigma }\frac{1}{\sigma _i^2}\right)}.$$
(29)
The sign reversal due to the field inversion was accounted for in the raw data, before the data were fit. Thus, the runs are combined using conventional (i.e., additive) averaging. The final sidereal amplitude $`A`$ was calculated by adding the mean cosine and sine amplitudes in quadrature, $`A=\sqrt{\overline{\delta \nu }_{Z,\alpha }^2+\overline{\delta \nu }_{Z,\beta }^2}`$. We measure a sidereal variation of the $`F=1`$, $`\mathrm{\Delta }m_F=\pm 1`$ Zeeman frequency of hydrogen of $`A=0.49\pm 0.34`$ mHz.

We note that since we are measuring an amplitude, and therefore a strictly positive quantity, this result is consistent with no sidereal variation at the 1-sigma level: in the case where $`\overline{\delta \nu }_{Z,\alpha }`$ and $`\overline{\delta \nu }_{Z,\beta }`$ have zero mean value and the same variance $`\sigma `$, the probability distribution for $`A`$ takes the form $`P(A)=A\sigma ^2\mathrm{exp}(A^2/2\sigma ^2)`$, which has the most probable value occurring at $`A=\sigma `$.

## V Error analysis

In addition to our automated acquisition of Zeeman frequency data, we continuosly monitored the maser’s external environment. At every ten second step, in addition to applied frequency and maser beat period, we recorded room temperature, maser cabinet temperature, solenoid current, maser power, ambient magnetic field, and active Helmholtz coil current (see Sec. V A).

### A Magnetic systematics

The $`F=1`$, $`m_F=\pm 1`$ Zeeman frequency depends to first-order on the z-component of the magnetic field in the storage bulb. Thus, all external field fluctuations must be sufficiently screened to enable a sensitivity to shifts from Lorentz and CPT symmetry violation. The maser cavity and bulb are therefore surrounded by a set of four nested magnetic shields that reduce the ambient field by a factor of about 32,000. We measure unshielded fluctuations in the ambient field of about 3 mG (peak-peak) during the day, and even when shielded, these add significant noise to a single Zeeman scan, as illustrated in Fig. 11(a). Furthermore, the amplitude of the field fluctuations is significantly reduced late at night, which could generate a diurnal systematic effect in our data.

To reduce the effect of fluctuations in the ambient magnetic field, we installed an active feedback system (see Fig. 10) consisting of two pairs of large Helmholtz coils (2.4 m diameter). The first pair of coils (50 turns) produced a uniform field that cancelled most of the z-component of the ambient field, leaving a residual field of around 5 mG. A magnetometer probe that sensed the residual ambient field was placed partially inside the maser’s magnetic shields near the maser cavity. This probe had a sensitivity of $`s`$ = 1.7 mG/V. Due to its location partially inside the magnetic shields, the probe was screened by a factor of about six from external fields, reducing the sensitivity to $`s^{}`$ = 0.3 mG/V, and producing a differential screening of 5300 between the magnetometer probe and the atoms. The magnetometer output was passed into a PID servo (Linear Research model LR-130), which contained a proportional stage (gain $`G`$ = 33), an integral stage (time constant Ti = 0.1 s) and a derivative stage (time constant = 0.01 s). The correction voltage was applied to the second pair of Helmholtz coils (3 turns) which produced a uniform field ($`p`$ = 14 mG/V) along the z-axis to nullify the residual field and actively counter any field fluctuations. Neglecting the small effect of the derivative stage, the overall time constant of this system was given by $`\tau =T_i(1+s^{}/pG)`$ 0.1 s, about 100 times shorter than the averaging time of our maser frequency shift measurements (10 s).

With this system we were able to further reduce ambient field fluctuations at the magnetometer by a factor of 3,000. The resulting unshielded fluctuations were less than 1 $`\mu `$G peak-peak. The field recorded by the partially screened magnetometer probe is shown in Fig. 12. The noise on a single Zeeman scan was reduced below our Zeeman frequency resolution, as shown in Fig. 11(b). During our Lorentz symmetry test, we monitored the field at the magnetometer probe and placed a bound of $``$ 5 nG on the sidereal component of the variation. This corresponds to a shift of less than 0.2 $`\mu `$Hz on the hydrogen Zeeman frequency, three orders of magnitude smaller than the sidereal Zeeman frequency bound measured.

The magnetometer  used in the feedback loop was a fluxgate magnetometer probe (RFL industries Model 101) which consisted of two parallel high-permeability magnetic cores each surrounded by an excitation coil (the excitation coils were wound in the opposite sense of each other). A separate pickup coil was wound around the pair of cores. An AC current (about 2.5 kHz) in the excitation coils drove the cores into saturation, and, in the presence of any slowly varying external magnetic field oriented along the magnetic cores’ axes, an EMF was generated in the pickup coil at the second and higher harmonics of the excitation frequency. The magnitude of the time-averaged EMF was proportional to the external field. The probe had a sensitivity of approximately 1 nG.

Any Lorentz violating spin-orientation dependence of the energy of the electrons in the magnetic cores would induce a sidereal variation in the cores’ magnetization and could generate, or mask, a sidereal variation in the hydrogen Zeeman frequency through the feedback circuit. However, based on the latest bound on electron Lorentz violation  (10-29 GeV), the Lorentz violating shift would be less than 10-11 G, far below the level of residual ambient field fluctuations. Also, the additional shielding factor of 5300 between the probe and the atoms further reduced the effect of any Lorentz violating shift in the probe electrons’ energies.

With the ambient field kept nearly constant near zero, the Zeeman frequency was set by the magnetic field generated by the solenoid, and hence by the solenoid current. We monitored solenoid current fluctuations by measuring the voltage across the current-setting 5 k$`\mathrm{\Omega }`$ resistor with a 5 1/2 digit multimeter (Fluke model 8840A/AF). By measuring the Zeeman frequency shift caused by large current changes, we found a dependence of around 10 mHz/nA. When acquiring Lorentz symmetry test data, we measured long term drifts in the current of about 5 nA (see Fig. 13), significant enough to produce detectable shifts in the Zeeman frequency. Thus, we subtracted these directly from the Zeeman data. We measured a sidereal variation of 25 $`\pm `$ 10 pA on the solenoid current, corresponding to a sidereal variation of 0.16 $`\pm `$ 0.08 mHz on the Zeeman frequency correction. This systematic uncertainty in the Zeeman frequency was included in the net error analysis, as described in Sec. V C.

### B Other systematics

The maser resided in a closed, temperature stabilized room where the temperature oscillated with a peak-peak amplitude of slightly less than 0.5˚C with a period of around 15 minutes. The maser was contained in an insulated and thermally controlled cabinet, which provided a factor of five to ten shielding from the room, and reduced the fluctuations to less than 0.1˚C peak-peak, as shown in Fig. 14. By making large changes in the maser cabinet temperature and measuring the effect on the Zeeman frequency, we found a temperature coefficient of about 200 mHz/˚C. We believe this frequency shift was due mainly to the resistors which set the solenoid current, which had 100 ppm/˚C temperature coefficients. We monitored the cabinet temperature and placed a bound on the sidereal component of the temperature fluctuations at 0.5 mK, which would produce a systematic sidereal variation of 100 $`\mu `$Hz on the Zeeman frequency, about a factor of 3 smaller than the measured limit on sidereal variation in Zeeman frequency.

As mentioned in Sec. III C, spin-exchange effects induce a small offset of the Zeeman frequency given by Eqn. 24 from the actual Zeeman frequency . This would imply that fluctuations in the input atomic flux (and therefore the maser power) could cause fluctuations in the Zeeman frequency measurement. We measured a limit on the shift of the Zeeman frequency due to large changes in average maser power at less than 0.8 mHz/fW. (Expected shifts from spin-exchange are ten times smaller than this level (Sec. III C). We believe the measured limit is related to heating of the maser as the flux is increased). During long-term operation, the average maser power drifted approximately 1 fW/day (see Fig. 15). The sidereal component of the variations of the maser power were less than 0.05 fW, implying a variation in the Zeeman frequency of less than 40 $`\mu `$Hz, an order of magnitude smaller than our experimental bound for sidereal Zeeman frequency variation.

### C Final result

We measured systematic errors in sidereal Zeeman frequency variation (as described in Secs. V A and V B) due to ambient magnetic field (0.2 $`\mu `$Hz), solenoid field (80 $`\mu `$Hz), maser cabinet temperature (100 $`\mu `$Hz), and hydrogen density induced spin-exchange shifts (40 $`\mu `$Hz). Combining these errors in quadrature with the 0.34 mHz statistical uncertainty in Zeeman frequency variation, we find a sidereal variation of the $`F=1`$, $`\mathrm{\Delta }m_F=\pm 1`$ Zeeman frequency in hydrogen of 0.44 $`\pm `$ 0.37 mHz at the 1-$`\sigma `$ level. This 0.37 mHz bound corresponds to 1.5 $`\times `$ 10-27 GeV in energy units.

## VI Discussion

### A Transformation to fixed frame

Our experimental bound of 0.37 mHz on sidereal variation of the hydrogen Zeeman frequency may be interpreted in terms of Eqn. 5 as a bound on vector and tensor components of the standard model extension. To make meaningful comparisons to other experiments, we transform our result into a fixed reference frame. Following the construction in reference , we label the fixed frame with coordinates (X,Y,Z) and the laboratory frame with coordinates (x,y,z), as shown in Fig. 16. We select the earth’s rotation axis as the fixed Z axis, (declination = 90 degrees). We then define fixed X as declination = right ascension = 0 degrees, and fixed Y as declination = 0 degrees, right ascension = 90 degrees. With this convention, the X and Y axes lie in the plane of the earth’s equator. Note that the $`\alpha `$, $`\beta `$ axes of Sec. IV E, also in the earth’s equatorial plane, are rotated about the earth’s rotation axis from the X,Y axes by an angle equivalent to the right ascension of 71 7’ longitude at 00:00 of November 19, 1999.

For our experiment, the quantization axis (which we denote z) was vertical in the lab frame, making an angle $`\chi `$ 48 degrees relative to Z, accounted for by rotating the entire (x,y,z) system by $`\chi `$ about Y. The lab frame (x,y,z) rotates about Z by an angle $`\mathrm{\Omega }t`$, where $`\mathrm{\Omega }`$ is the frequency of the earth’s (sidereal) rotation.

These two coordinate systems are related through the transformation

$$\left(\begin{array}{c}t\\ x\\ y\\ z\end{array}\right)=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& \mathrm{cos}\chi \mathrm{cos}\mathrm{\Omega }t& \mathrm{cos}\chi \mathrm{sin}\mathrm{\Omega }t& \mathrm{sin}\chi \\ 0& \mathrm{sin}\mathrm{\Omega }t& \mathrm{cos}\mathrm{\Omega }t& 0\\ 0& \mathrm{sin}\chi \mathrm{cos}\mathrm{\Omega }t& \mathrm{sin}\chi \mathrm{sin}\mathrm{\Omega }t& \mathrm{cos}\chi \end{array}\right)\left(\begin{array}{c}0\\ X\\ Y\\ Z\end{array}\right)=𝐓\left(\begin{array}{c}0\\ X\\ Y\\ Z\end{array}\right).$$
(30)
Then, vectors transform as $`\stackrel{}{b}_{lab}=𝐓\stackrel{}{b}_{fixed}`$, while tensors transform as $`𝐝_{lab}=𝐓𝐝_{fixed}𝐓^1`$.

As shown in equation (5), our signal depends on the following combination of terms (for both electron and proton):

$$\stackrel{~}{b}_3=b_3md_{30}H_{12}.$$
(31)
Transforming these to the fixed frame, we see

$`b_3`$ $`=`$ $`b_Z\mathrm{cos}\chi +b_X\mathrm{sin}\chi \mathrm{cos}\mathrm{\Omega }t+b_Y\mathrm{sin}\chi \mathrm{sin}\mathrm{\Omega }t,`$ (32)
$`d_{30}`$ $`=`$ $`d_{Z0}\mathrm{cos}\chi +d_{X0}\mathrm{sin}\chi \mathrm{cos}\mathrm{\Omega }t+d_{Y0}\mathrm{sin}\chi \mathrm{sin}\mathrm{\Omega }t,`$ (33)
$`H_{12}`$ $`=`$ $`H_{XY}\mathrm{cos}\chi +H_{YZ}\mathrm{sin}\chi \mathrm{cos}\mathrm{\Omega }t+H_{ZX}\mathrm{sin}\chi \mathrm{sin}\mathrm{\Omega }t,`$ (34)
so our observable is given by

$`\stackrel{~}{b}_3`$ $`=`$ $`\left(b_Zmd_{Z0}H_{XY}\right)\mathrm{cos}\chi `$ (35)
$`+`$ $`\left(b_Ymd_{Y0}H_{ZX}\right)\mathrm{sin}\chi \mathrm{sin}\mathrm{\Omega }t`$ (36)
$`+`$ $`\left(b_Xmd_{X0}H_{YZ}\right)\mathrm{sin}\chi \mathrm{cos}\mathrm{\Omega }t.`$ (37)
The first term on the right is a constant offset, not bounded by our experiment. The second and third terms each vary at the sidereal frequency. Combining Eqn. 35 (for both e- and p) with Eqn. 5, we see

$`|\mathrm{\Delta }\nu _Z|^2`$ $`=`$ $`\left[\left(b_Y^em_ed_{Y0}^eH_{ZX}^e\right)+\left(b_Y^pm_pd_{Y0}^pH_{ZX}^p\right)\right]^2{\displaystyle \frac{\mathrm{sin}^2\chi }{h^2}}`$ (38)
$`+`$ $`\left[\left(b_X^em_ed_{X0}^eH_{YZ}^e\right)+\left(b_X^pm_pd_{X0}^pH_{YZ}^p\right)\right]^2{\displaystyle \frac{\mathrm{sin}^2\chi }{h^2}}.`$ (39)
Inserting $`\chi `$ = 48 degrees, we obtain the final result

$$\sqrt{\left(\stackrel{~}{b}_X^e+\stackrel{~}{b}_X^p\right)^2+\left(\stackrel{~}{b}_Y^e+\stackrel{~}{b}_Y^p\right)^2}=\left(3\pm 2\right)\times 10^{27}GeV.$$
(40)
Our 1-sigma bound on Lorentz and CPT violation of the proton and electron is therefore 2 $`\times `$ 10-27 GeV.

### B Comparison to previous experiments

We compare our result with other recent tests of Lorentz and CPT symmetry in Table III. Although our bounds are numerically similar to the those from the 199Hg/133Cs experiment, the simplicity of the hydrogen atom allows us to place bounds directly on the electron and proton; uncertainties in nuclear structure models do not complicate the interpretation of our result. The recent limit set by the torsion pendulum experiment of Adelberger et. al.  on electron Lorentz and CPT violation casts our result as a clean bound on Lorentz and CPT violation of the proton.

### C Future work

To make a more sensitive measure of the sidereal variation of the Zeeman frequency in a hydrogen maser, it will be important to clearly identify and reduce the magnitude of the long term drifts of the Zeeman frequency. Possible sources of these drifts are magnetic fields near the maser bulb caused by stray currents in heaters or power supplies in the inner regions of the maser. Also, the scatter of the Zeeman data points, believed to be due mainly to residual thermal fluctuations, should be reduced. Both of these objectives could be accomplished by carefully rebuilding a hydrogen maser, with better engineered power and temperature control systems.

## VII Acknowledgments

We gratefully acknowledge the encouragement of Alan Kostelecký. Financial support was provided by NASA grant NAG8-1434 and ONR grant N00014-99-1-0501. M.A.H. thanks NASA for fellowship support under the Graduate Student Researchers Program.


text:
# A Transition to the Soft State in GRS 1758–258

## 1 Introduction

GRS 1758–258 is one of only three black-hole candidates (BHCs) usually near its maximum luminosity and usually in the hard state (the other two are  Cyg X–1 and  1E 1740.7–2942). It displays shot-noise flickering and quasi-periodic oscillations with characteristic timescales on the order of 1 s (Smith et al., 1997; Lin et al., 2000). It shows a core source and lobes in the radio (Rodriguez et al., 1992), which caused it to be one of the first sources referred to as a “microquasar”.

The spectrum of  GRS 1758–258 is generally a power law (with an exponential cutoff beginning around 100 keV), occasionally with a weak thermal component (Mereghetti et al., 1994; Heindl & Smith, 1998). Weekly observations with the Proportional Counter Array (PCA) (Jahoda et al., 1996) on the Rossi X-ray Timing Explorer (RXTE) show that the power-law index can vary from about 1.5 to 2.5, that the hard x-ray luminosity is always dominated by the power law, and that the variability before 2001 (Figure 1) was not extreme (Main et al., 1999; Smith et al., 2001a). Earlier observations by GRANAT/Sigma at 40-150 keV showed the source below detection level in Fall 1991 and Spring 1992 (Gilfanov et al., 1993), and therefore more than an order of magnitude fainter than its typical GRANAT and RXTE luminosity. The Burst and Transient Source Experiment (BATSE) on the Compton Gamma-Ray Observatory also saw the 20-100 keV flux from  GRS 1758–258 vary dramatically (Zhang et al., 1997c).

GRS 1758–258’s spectral variations are not simultaneous with variations in its luminosity, a newly-discovered characteristic which it shares with  1E 1740.7–2942 (Main et al., 1999; Smith et al., 2001a). In  Cyg X–1, the usual prototype of a BHC system, there is a clearly defined progression from the hard to soft state: the 2-10 keV flux increases dramatically as the spectrum softens, the total x-ray luminosity increasing only slightly (Zhang et al., 1997a; Gierliński et al., 1999). The softening takes two forms: the index of the power law softens while at the same time the thermal component brightens and goes to higher temperature. In the transition back to the hard state, the spectrum and luminosity also vary simultaneously.

GRS 1758–258 and  1E 1740.7–2942, in contrast, tend to be softest not when the 2-10 keV flux is highest, but rather when the derivative of that count rate is most negative (Smith et al., 2001a). We interpreted this behavior in terms of a model by Chakrabarti & Titarchuk (1995) in which BHC accretion operates in two decoupled flows, a classical thin disk which produces thermal emission and a nearly spherical, sub-Keplerian flow (“halo” hereafter) which scatters the thermal photons into the cut-off power law. If mass transfer from the companion suddenly dropped, altering both flows at large radii, the change would propagate at the free-fall timescale in the halo, but at the much slower viscous timescale in the thin disk. Therefore, there would temporarily be an imbalance: the same number of soft disk photons, but fewer hot halo electrons to scatter them. Fewer scatters per photon softens the power law and, in addition, the thinner halo would be more easily cooled by the disk photons. In Smith et al. (2001a) we dubbed this a “dynamical” soft state and hypothesized that it was responsible for almost all the softening events observed in  GRS 1758–258 and  1E 1740.7–2942, except for a brief period in Fall of 2000 when  GRS 1758–258 was at its brightest and seemed on the verge of entering a true (“static”) soft state in which the halo self-cools and collapses.

In this Letter we report on the recent transition of  GRS 1758–258 to the soft state, and interpret the data in terms of the dynamical soft state model.

## 2 Observations

PCA snapshots of 1500 s were taken monthly of  GRS 1758–258 in 1996, weekly through 2000, and are being taken twice weekly from March 2001. There are no pointings for a period from November to January each year due to a Sun-angle constraint on pointing. Figure 1 shows the count rate as a function of time in two energy bands: 2.5-4.0 keV and 10.0-25.0 keV. Instrumental background and background due to Galactic Plane diffuse emission (Smith et al., 1997; Main et al., 1999) have been subtracted. The source is kept off-axis to prevent contamination from GX 5-1 nearby. All results below use data from layer 1 of the PCA only. The data were analyzed with version 5.0.4 of HEAsoft.

There was a precipitous drop above 10 keV in late February of 2001 (Figure 1), but not in the soft band (Smith et al., 2001b, c). On March 12-13, 31 ksec of public observations were made with RXTE in three separate intervals. Since the new soft state would have been undetectable to Sigma, it is quite possible that this was the state that instrument observed in 1991-1992.

## 3 Results

Figure 2 shows the PCA spectrum of the combined March 12-13 data, including the residuals from a fit to a power law plus a blackbody. Because the PCA fits begin at 3 keV, the fitted hydrogen absorption column and blackbody intensity are strongly correlated and can’t both be determined well. We therefore fixed the column at the ASCA value (Mereghetti et al., 1997): 1.5 $`\times 10^{22}`$ atoms cm-2. The fitted blackbody temperature was (0.395 $`\pm `$ 0.006) keV and the power-law photon index was (2.89 $`\pm `$ 0.12). Using a multicolor disk blackbody, the maximum temperature was (0.464 $`\pm `$ 0.007) keV and the power law index (2.75 $`\pm `$ 0.12). The blackbody luminosity ($`3.0\times 10^9`$ ergs cm-2s-1, or $`2.6\times 10^{37}`$ ergs s-1 at 8.5 kpc) is the highest seen with RXTE , although only about 35% higher than in some spectra from near the peak in the top panel of Figure 1 in late 2000.

As is apparent from the residuals in Figure 2, the addition of an iron K-shell line improves the fit. Fixing the energy at 6.7 keV and the width at 0.1 keV (narrower than the PCA resolution), the derived equivalent width is 300 eV for both the blackbody and disk blackbody models. The interpretation of a weak, instrumentally broad line near the crossover point of two continuum components is problematic, and we do not claim this as a certain detection. No significant improvement in $`\chi ^2`$ resulted from letting the energy or width of the line vary.

Figure 3 shows the fits to the monitoring observations in early 2001, before and after the state transition. The first observation of 2001 (January 23, labeled “A” in Figure 3) displays both a bright power law and a significant blackbody peak (this is also occasionally seen at other periods from 1996-2000). The soft excess from 1.0-2.4 keV measured by Mereghetti et al. (1994) with ROSAT in 1994 was 2.8$`\times 10^{10}`$ ergs cm-2s-1, about half the value for this spectrum (5.7$`\times 10^{10}`$ ergs cm-2s-1 from 1.0-2.4 keV, or 7.0$`\times 10^{10}`$ ergs cm-2s-1 bolometric). The remaining spectra before the transition have no detectable blackbody and a harder power-law index. The first spectrum after the transition (labeled “C”) has somewhat more of a residual power-law tail than the spectra after it.

Figure 4a shows a combined lightcurve around the time of the transition, using data from two sources: the monitoring campaign discussed above, and a second RXTE campaign of twice-weekly scans of the whole Galactic-bulge region (Markwardt et al., 1999, 2001). For each scan, the data are fitted to a model of point sources and diffuse emission. One scan point on February 20 implies that the source may have dropped to a low level in less than 3.9 dy and then rebounded completely less than 11.3 hr later. Generally speaking the pointed observations and bulge-scan intensities agree to within $`2\sigma `$. However, the proximity of  GRS 1758–258 to the bright source GX 5–1 may potentially introduce unmodelled systematic errors into the scan data.

Figure 4b shows the evolution of the root mean square (rms) variability beyond Poisson noise. As is typical for BHC soft states (e.g. Maejima et al., 1984), the variability drops dramatically. The last data point represents all the deep-pointing data from March 12-13. At this time the rms variability (0.004-16 Hz, 2-10 keV) was consistent with zero, with a 3$`\sigma `$ upper limit of 2.5%.

## 4 Discussion

There are similarities and differences between this hard-to-soft transition and transitions in other BHCs. One of the best-monitored transitions was in  Cyg X–1 in 1996. The 1.3-200 keV luminosity as measured by the RXTE All-Sky Monitor (ASM) and BATSE, corrected for interstellar absorption, jumped upwards by about 35% (Zhang et al., 1997a). To compare, we extrapolated our spectral fits to this energy range and removed the interstellar absorption using the fixed column density discussed above.  GRS 1758–258 and other BHCs have an exponential cutoff in the hard state beyond the PCA range, and we used cutoff parameters from a deep RXTE pointing to  GRS 1758–258 in 1996 (Heindl & Smith, 1998). We find that the extrapolated 1.3-200 keV luminosity of  GRS 1758–258 dropped in the hard-to-soft transition by about the same amount that it rose in  Cyg X–1 in 1996.

That soft state in  Cyg X–1, however, was not typical of BHC soft states in general. The power-law flux at 25 keV dropped only by a factor of 3 from the hard state (Zhang et al., 1997b) compared to the factor of 50 apparent in Figure 3, and there was still 20% rms of fast variability (Cui et al., 1997). This led Belloni et al. (1996) to conclude that this “soft” state should have been called an “intermediate” state. If the hard/intermediate/soft progression is a monotonic increase in accretion rate as generally thought, one would expect that the more “complete” transition in  GRS 1758–258 would result in an even larger increase in luminosity, not a decrease.

GX 339–4 entered the soft state in 1990 and 1998 with similar spectral and timing characteristics to the new soft state in  GRS 1758–258. Between April and August of 1990, the power-law tail in  GX 339–4 dropped by at least an order of magnitude at 10 keV and the rms variability in the soft state was only (6.1 $`\pm `$ 2.7)% (Grebenev et al., 1991). For the 1998 transition, Belloni et al. (1999) compared the 2.5-20 keV unabsorbed flux in the soft state to data taken in the hard state several months previously. They found that the total flux in this band increased by about a factor of two. We can evaluate our PCA data in this band without extrapolation. We find again that, contrary to the behavior of the other BHC,  GRS 1758–258 saw a decrease in absorption-corrected flux in the energy band chosen for comparison, in this case a drop of a factor of 3.3 from 2.5-20 keV between February 21 and March 2.

BATSE was the first to observe a soft state in  1E 1740.7–2942, with a power law index of -2.6 from 20-100 keV (Zhang et al., 1997c). The soft spectrum crossed the usual harder spectrum at 20 keV, so again there would have been an increase rather than a drop in the power-law component over the 3-25 keV band.

While the transitions in  Cyg X–1,  GX 339–4, and other BHCs are probably indeed caused by an increase in accretion causing self-cooling and partial collapse of the halo (static soft state), we suggest that this transition in  GRS 1758–258 is instead caused by the sudden shutoff of all, or nearly all, of the mass transfer from the compact object’s companion. This is the extreme instance of the dynamical soft state described in Smith et al. (2001a). Given the hypothesis (Chakrabarti & Titarchuk, 1995) that the disk and halo flows are independently fed directly from the companion, the halo would then vanish immediately leaving the disk to decay away on a viscous timescale on the order of a month or more (Main et al., 1999; Smith et al., 2001a). We wouldn’t expect  Cyg X–1 to show a dynamical soft state, because as a wind accretor its disk is expected to be small, and therefore have a short viscous time (Smith et al., 2001a).

It is possible that, beneath the halo, the disk emission was more or less constant throughout early 2001, and when it wasn’t seen, the halo was optically thick. In the first spectrum of the year (“A” in Figure 3), with a blackbody luminosity about a factor of 4 below the naked disk, the optical depth of the halo would then have been near unity. The changes in the slope of the power law support this interpretation: the slope is hardest when the halo is assumed to be thickest, and repeated scatterings should produce a harder power law (Shapiro et al., 1976). If the apparent drop on February 20 (Figure 4a) is genuine, then the time for the halo to re-establish itself would be less than 11.3 hr, which is still longer than the free-fall time from the outer disk.

Viewing the disk and halo flows as independent, and assuming that the disk flow did not change as the halo disappeared, the loss of fast variability in the transition implies that the flickering commonly observed in the hard state of BHCs is either intrinsic to the halo flow alone, or else due to an interaction at the boundary between the dense, radially slow-moving disk and the thin, radially fast-moving halo.

When the transition to the soft state occurred, we predicted (Smith et al., 2001c), based on the dynamical picture and the disk timescales we previously derived (Main et al., 1999; Smith et al., 2001a), that mass transfer from the companion had stopped and the soft-state emission would decay away with a timescale of a few weeks. Figure 5 shows the flux (uncorrected for absorption) from 3-25 keV as a function of time starting with the first soft-state spectrum (February 27) to the most recent data available (April 15). The shape is neither exponential or linear, but somewhere in between. The best fit exponential gives a (27.9 $`\pm `$ 1.7) dy time constant. The latest fluxes approach the limit of the systematic uncertainty in the Galactic diffuse emission.

If there is still a small mass transfer from the companion, the spectrum may eventually make a transition to the pure, very faint hard power law seen in  GX 339–4 in its lowest luminosity (“off”) state (Maejima et al., 1984; Ilovaisky et al., 1986; Mendez & van der Klis, 1997; Asai et al., 1998; Kong et al., 2000). Observations of  GX 339–4 were never taken immediately following the drop from higher fluxes, and it may show a dynamical soft state at those times. Eventually, our understanding of the state change in  GRS 1758–258 may apply to the state changes in the soft x-ray transients (the outbursts of the more numerous BHCs without persistent emission), which can show hysteresis with a hard-to-soft transition near peak flux followed by a much later soft-to-hard transition at much lower flux (Miyamoto et al., 1995).