From 02187d60a03d2f7bed87813fe4c15646c93b5eb5 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Sat, 11 Apr 2020 20:54:10 -0400
Subject: Added Appendix with calculation of modulus with higher-order
 interaction, made figure bigger.

---
 main.tex | 287 ++++++++++++++++++++++++++++++++++++++++++++++++++++++---------
 1 file changed, 248 insertions(+), 39 deletions(-)

(limited to 'main.tex')

diff --git a/main.tex b/main.tex
index 7e934bb..b815c79 100644
--- a/main.tex
+++ b/main.tex
@@ -1,5 +1,5 @@
 
-\documentclass[aps,prb,reprint,longbibliography,floatfix]{revtex4-1}
+\documentclass[aps,prb,reprint,longbibliography,floatfix,fleqn]{revtex4-1}
 \usepackage[utf8]{inputenc}
 \usepackage{amsmath,graphicx,upgreek,amssymb,xcolor}
 \usepackage[colorlinks=true,urlcolor=purple,citecolor=purple,filecolor=purple,linkcolor=purple]{hyperref}
@@ -111,9 +111,10 @@ broken symmetry remains unknown. This state, known as \emph{hidden order}
 lower temperatures.  At sufficiently large hydrostatic pressures, both
 superconductivity and \ho\ give way to local moment antiferromagnetism
 (\afm).\cite{Hassinger_2008}  Modern theories~\cite{Kambe_2018, Haule_2009,
-  Kusunose_2011, Kung_2015, Cricchio_2009, Ohkawa_1999, Santini_1994,
-Kiss_2005, Harima_2010, Thalmeier_2011, Tonegawa_2012, Rau_2012, Riggs_2015,
-Hoshino_2013, Ikeda_1998, Chandra_2013a, 1902.06588v2, Ikeda_2012} propose
+Kusunose_2011_On, Kung_2015, Cricchio_2009, Ohkawa_1999, Santini_1994,
+Kiss_2005, Harima_2010, Thalmeier_2011, Tonegawa_2012_Cyclotron,
+Rau_2012_Hidden, Riggs_2015_Evidence, Hoshino_2013_Resolution,
+Ikeda_1998_Theory, Chandra_2013_Hastatic, 1902.06588v2, Ikeda_2012} propose
 associating any of a variety of broken symmetries with \ho. This work analyzes
 a family of phenomenological models with order parameters of general symmetry
 that couple linearly to strain. Of these, only one is compatible with two
@@ -128,9 +129,10 @@ examine the thermodynamic discontinuities in the elastic moduli at
 $T_{\text{\ho}}$.\cite{1903.00552v1} The observation of discontinues only in
 compressional, or $\Aog$, elastic moduli requires that the point-group
 representation of \ho\ be one-dimensional. This rules out many order parameter
-candidates~\cite{Thalmeier_2011, Tonegawa_2012, Rau_2012, Riggs_2015,
-Hoshino_2013, Ikeda_2012, Chandra_2013b} in a model-independent way, but
-doesn't differentiate between those that remain. 
+candidates~\cite{Thalmeier_2011, Tonegawa_2012_Cyclotron, Rau_2012_Hidden,
+Riggs_2015_Evidence, Hoshino_2013_Resolution, Ikeda_2012, Chandra_2013_Origin}
+in a model-independent way, but doesn't differentiate between those that
+remain. 
 
 Recent x-ray experiments discovered rotational symmetry breaking in \urusi\
 under pressure.\cite{Choi_2018} Above 0.13--0.5 $\GPa$ (depending on
@@ -462,33 +464,34 @@ signatures of a continuous transition by locating thermodynamic singularities
 at nonzero $q=q_*$.  The remaining clue at $q=0$ is a particular kink in the
 corresponding modulus.
 
-\begin{figure}[htpb]
+\section{Comparison to experiment}
+
+\begin{figure*}[htpb]
   \centering
-  \includegraphics[width=\columnwidth]{fig-stiffnesses}
+  \includegraphics{fig-stiffnesses}
   \caption{
-   \Rus\ measurements of the elastic moduli of \urusi\ at ambient pressure as a
-   function of temperature from recent experiments\cite{1903.00552v1} (blue,
-   solid) alongside fits to theory (magenta, dashed). The solid yellow region
-   shows the location of the \ho\ phase. (a) $\Btg$ modulus data and a fit to
-   the standard form.\cite{Varshni_1970} (b) $\Bog$ modulus data and a fit to
-   \eqref{eq:static_modulus}. The fit gives
-   $C^0_\Bog\simeq\big[71-(0.010\,\K^{-1})T\big]\,\GPa$, $D_\perp
-   q_*^4/b^2\simeq0.16\,\GPa^{-1}$, and
-   $a/b^2\simeq6.1\times10^{-4}\,\GPa^{-1}\,\K^{-1}$. Addition of a quadratic
-   term in $C^0_\Bog$ was here not needed for the fit.\cite{Varshni_1970} (c)
-   $\Bog$ modulus data and the fit of the \emph{bare} $\Bog$ modulus. (d)
-   $\Bog$ modulus data and the fit transformed by
-   $[C^0_\Bog(C^0_\Bog/C_\Bog-1)]]^{-1}$, which is predicted from
-   \eqref{eq:static_modulus} to equal $D_\perp q_*^4/b^2+a/b^2|T-T_c|$, e.g.,
-   an absolute value function. The failure of the Ginzburg--Landau prediction
-   below the transition is expected on the grounds that the \op\ is too large
-   for the free energy expansion to be valid by the time the Ginzburg
-   temperature is reached.
+    \Rus\ measurements of the elastic moduli of \urusi\ at ambient pressure as a
+    function of temperature from recent experiments\cite{1903.00552v1} (blue,
+    solid) alongside fits to theory (magenta, dashed). The solid yellow region
+    shows the location of the \ho\ phase. (a) $\Btg$ modulus data and a fit to
+    the standard form.\cite{Varshni_1970} (b) $\Bog$ modulus data and a fit to
+    \eqref{eq:static_modulus}. The fit gives
+    $C^0_\Bog\simeq\big[71-(0.010\,\K^{-1})T\big]\,\GPa$, $D_\perp
+    q_*^4/b^2\simeq0.16\,\GPa^{-1}$, and
+    $a/b^2\simeq6.1\times10^{-4}\,\GPa^{-1}\,\K^{-1}$. Addition of a quadratic
+    term in $C^0_\Bog$ was here not needed for the fit.\cite{Varshni_1970} (c)
+    $\Bog$ modulus data and the fit of the \emph{bare} $\Bog$ modulus. (d)
+    $\Bog$ modulus data and the fit transformed by
+    $[C^0_\Bog(C^0_\Bog/C_\Bog-1)]]^{-1}$, which is predicted from
+    \eqref{eq:static_modulus} to equal $D_\perp q_*^4/b^2+a/b^2|T-T_c|$, e.g.,
+    an absolute value function. The failure of the Ginzburg--Landau prediction
+    below the transition is expected on the grounds that the \op\ is too large
+    for the free energy expansion to be valid by the time the Ginzburg
+    temperature is reached.
   }
   \label{fig:data}
-\end{figure}
+\end{figure*}
 
-\section{Comparison to experiment}
 \Rus\ experiments~\cite{1903.00552v1} yield the individual elastic moduli
 broken into irreps; data for the $\Bog$ and $\Btg$ components defined in
 \eqref{eq:strain-components} are shown in Figures \ref{fig:data}(a--b).  The
@@ -497,18 +500,27 @@ the presence of the transition, exhibiting the expected linear stiffening upon
 cooling from room temperature, with a low-temperature cutoff at some fraction
 of the Debye temperature.\cite{Varshni_1970} The $\Bog$ modulus
 Fig.~\ref{fig:data}(b) has a dramatic response, softening over the course of
-roughly $100\,\K$ and then cusping at the \ho\ transition. While the
-low-temperature response is not as dramatic as the theory predicts, mean field
-theory---which is based on a small-$\eta$ expansion---will not work
-quantitatively far below the transition where $\eta$ has a large nonzero value
-and higher powers in the free energy become important. The data in the
+roughly $100\,\K$ and then cusping at the \ho\ transition. The data in the
 high-temperature phase can be fit to the theory \eqref{eq:static_modulus}, with
 a linear background modulus $C^0_\Bog$ and $\tilde r-\tilde r_c=a(T-T_c)$, and
-the result is shown in Figure \ref{fig:data}(b). The data and theory appear
-quantitatively consistent in the high temperature phase, suggesting that \ho\
-can be described as a $\Bog$-nematic phase that is modulated at finite $q$
-along the $c-$axis. The predicted softening appears over hundreds of Kelvin;
-Figures \ref{fig:data}(c--d) show the background modulus $C_\Bog^0$ and the
+the result is shown in Figure \ref{fig:data}(b).
+
+The behavior of the modulus below the transition does not match
+\eqref{eq:static_modulus} well, but this is because of the truncation of the
+free energy expansion used above.  Higher order terms like $\eta^2\epsilon^2$
+contribute to the modulus starting at order $\eta_*^2$, and therefore while
+they do not affect the behavior above the transition, they change the behavior
+below it. To demonstrate this, in Appendix~\ref{sec:higher-order} we compute
+the modulus in a theory where the interaction free energy is truncated after
+fourth order with new term $\frac12g\eta^2\epsilon^2$. The thin solid black
+line in Fig.~\ref{fig:data} shows the fit of the \rus\ data to \eqref{eq:C0}
+and shows that high-order corrections can account for the low-temperature
+behavior. 
+
+The data and theory appear quantitatively consistent, suggesting that \ho\ can
+be described as a $\Bog$-nematic phase that is modulated at finite $q$ along
+the $c-$axis. The predicted softening appears over hundreds of Kelvin; Figures
+\ref{fig:data}(c--d) show the background modulus $C_\Bog^0$ and the
 \op--induced response isolated from each other. 
 
 We have seen that the mean-field theory of a $\Bog$ \op\ recreates the topology
@@ -624,6 +636,203 @@ such as ultrasound, that could further support or falsify this idea.
   Elena Hassinger. We thank Sayak Ghosh for \rus\ data.
 \end{acknowledgements}
 
+\appendix
+
+\section{Adding a higher-order interaction}
+\label{sec:higher-order}
+
+In this appendix, we compute the $\Bog$ modulus for a theory with a high-order
+interaction truncation to better match the low-temperature behavior.  Consider
+the free energy density $f=f_\e+f_\i+f_\op$ with
+\begin{equation}
+  \begin{aligned}
+    f_\e&=\frac12C_0\epsilon^2 \\
+    f_\i&=-b\epsilon\eta+\frac12g\epsilon^2\eta^2 \\
+    f_\op&=\frac12\big[r\eta^2+c_\parallel(\nabla_\parallel\eta)^2+c_\perp(\nabla_\perp\eta)^2+D(\nabla_\perp^2\eta)^2\big]+u\eta^4.
+  \end{aligned}
+  \label{eq:new_free_energy}
+\end{equation}
+The mean-field stain conditioned on the order parameter is found from
+\begin{equation}
+  \begin{aligned}
+    0
+    &=\frac{\delta F[\eta,\epsilon]}{\delta\epsilon(x)}\bigg|_{\epsilon=\epsilon_\star[\eta]} \\
+    &=C_0\epsilon_\star[\eta](x)-b\eta(x)+g\epsilon_\star[\eta](x)\eta(x)^2,
+  \end{aligned}
+\end{equation}
+which yields
+\begin{equation}
+  \epsilon_\star[\eta](x)=\frac{b\eta(x)}{C_0+g\eta(x)^2}.
+  \label{eq:epsilon_star}
+\end{equation}
+Upon substitution into \eqref{eq:new_free_energy} and expanded to fourth order
+in $\eta$, $F[\eta,\epsilon_\star[\eta]]$ can be written in the form
+$F_\op[\eta]$ alone with $r\to\tilde r=r-b^2/C_0$ and $u\to\tilde
+u=u+b^2g/2C_0^2$. The phase diagram in $\eta$ follows as before with the
+shifted coefficients, and namely $\langle\eta(x)\rangle=\eta_*\cos(q_*x_3)$ for
+$\tilde r<c_\perp^2/4D=\tilde r_c$ with $q_*^2=-c_\perp/2D$ and
+\begin{equation}
+  \eta_*^2=\frac{c_\perp^2-4D\tilde r}{12D\tilde u}
+  =\frac{|\Delta\tilde r|}{3\tilde u}.
+\end{equation}
+We would like to calculate the $q$-dependant modulus
+\begin{equation}
+  C(q)
+  =\frac1V\int dx\,dx'\,C(x,x')e^{-iq(x-x')},
+\end{equation}
+where
+\begin{widetext}
+\begin{equation}
+  C(x,x')
+  =\frac{\delta^2F[\eta_\star[\epsilon],\epsilon]}{\delta\epsilon(x)\delta\epsilon(x')}\bigg|_{\epsilon=\langle\epsilon\rangle}
+    =\frac{\delta^2F_\e[\eta_\star[\epsilon],\epsilon]}{\delta\epsilon(x)\delta\epsilon(x')}+
+     \frac{\delta^2F_\i[\eta_\star[\epsilon],\epsilon]}{\delta\epsilon(x)\delta\epsilon(x')}+
+     \frac{\delta^2F_\op[\eta_\star[\epsilon],\epsilon]}{\delta\epsilon(x)\delta\epsilon(x')}
+     \bigg|_{\epsilon=\langle\epsilon\rangle}
+\end{equation}
+and $\eta_\star$ is the mean-field order parameter conditioned on the strain defined implicitly by
+\begin{equation}
+  0=\frac{\delta F[\eta,\epsilon]}{\delta\eta(x)}\bigg|_{\eta=\eta_\star[\epsilon]}
+  =-b\epsilon(x)+g\epsilon(x)^2\eta_\star[\epsilon](x)+\frac{\delta F_\op[\eta]}{\delta\eta(x)}\bigg|_{\eta=\eta_\star[\epsilon]}.
+  \label{eq:eta_star}
+\end{equation}
+We will work this out term by term. The elastic term is the most straightforward, giving
+\begin{equation}
+  \frac{\delta^2F_\e[\epsilon]}{\delta\epsilon(x)\delta\epsilon(x')}
+  =\frac12C_0\frac{\delta^2}{\delta\epsilon(x)\delta\epsilon(x')}\int dx''\,\epsilon(x'')^2
+  =C_0\delta(x-x').
+\end{equation}
+The interaction term gives
+\begin{equation}
+  \begin{aligned}
+    \frac{\delta^2F_\i[\eta_\star[\epsilon],\epsilon]}{\delta\epsilon(x)\delta\epsilon(x')}
+    &=-b\frac{\delta^2}{\delta\epsilon(x)\delta\epsilon(x')}\int dx''\,\epsilon(x'')\eta_\star[\epsilon](x'')
+      +\frac12g\frac{\delta^2}{\delta\epsilon(x)\delta\epsilon(x')}\int dx''\,\epsilon(x'')^2\eta_\star[\epsilon](x'')^2 \\
+    &=-b\frac{\delta\eta_\star[\epsilon](x')}{\delta\epsilon(x)}
+    -b\frac{\delta}{\delta\epsilon(x)}\int dx''\,\epsilon(x'')\frac{\delta\eta_\star[\epsilon](x'')}{\delta\epsilon(x')}
+    +g\frac{\delta}{\delta\epsilon(x)}\big[\epsilon(x')\eta_\star[\epsilon](x')^2\big] \\
+    &\qquad+g\frac{\delta}{\delta\epsilon(x)}\int dx''\,\epsilon(x'')^2\eta_\star[\epsilon](x'')\frac{\delta\eta_\star[\epsilon](x'')}{\delta\epsilon(x')} \\
+    &=-2(b-2g\epsilon(x)\eta_\star[\epsilon](x))\frac{\delta\eta_\star[\epsilon](x)}{\delta\epsilon(x')}-b\int dx''\,\epsilon(x'')\frac{\delta^2\eta_\star[\epsilon](x'')}{\delta\epsilon(x)\delta\epsilon(x')}
+    +g\eta_\star[\epsilon](x)^2\delta(x-x') \\
+    &\qquad+g\int dx''\,\epsilon(x'')^2\frac{\delta\eta_\star[\epsilon](x'')}{\delta\epsilon(x)}\frac{\delta\eta_\star[\epsilon](x'')}{\delta\epsilon(x')}
+    +g\int dx''\,\epsilon(x'')^2\eta_\star[\epsilon](x'')\frac{\delta^2\eta_\star[\epsilon](x'')}{\delta\epsilon(x)\delta\epsilon(x')}.
+  \end{aligned}
+\end{equation}
+The order parameter term relies on some other identities. First, \eqref{eq:eta_star} implies
+\begin{equation}
+  \frac{\delta F_\op[\eta]}{\delta\eta(x)}\bigg|_{\eta=\eta_\star[\epsilon]}
+    =b\epsilon(x)-g\epsilon(x)^2\eta_\star[\epsilon](x),
+  \label{eq:dFodeta}
+\end{equation}
+and therefore that the functional inverse $\eta_\star^{-1}[\eta]$ is
+\begin{equation}
+  \eta_\star^{-1}[\eta](x)=\frac b{2g\eta(x)}\Bigg(1-\sqrt{1-\frac{4g\eta(x)}{b^2}\frac{\delta F_\op[\eta]}{\delta\eta(x)}}\Bigg).
+\end{equation}
+The inverse function theorem further implies (with substitution of \eqref{eq:dFodeta} after the derivative is evaluated) that
+\begin{equation}
+  \bigg(\frac{\delta\eta_\star[\epsilon](x)}{\delta\epsilon(x')}\bigg)^{\{-1\}}
+  =\frac{\delta\eta_\star^{-1}[\eta](x)}{\delta\eta(x')}\bigg|_{\eta=\eta_\star[\epsilon]}
+  =\frac{g\epsilon(x)^2\delta(x-x')+\frac{\delta^2F_\op[\eta]}{\delta\eta(x)\delta\eta(x')}\big|_{\eta=\eta_\star[\epsilon]}}{b-2g\epsilon(x)\eta_\star[\epsilon](x)}
+\end{equation}
+and therefore that
+\begin{equation}
+  \frac{\delta^2F_\op[\eta]}{\delta\eta(x)\delta\eta(x')}\bigg|_{\eta=\eta_\star[\epsilon]}
+  =(b-2g\epsilon(x)\eta_\star[\epsilon](x))\bigg(\frac{\delta\eta_\star[\epsilon](x)}{\delta\epsilon(x')}\bigg)^{\{-1\}}
+  -g\epsilon(x)^2\delta(x-x').
+  \label{eq:d2Fodetadeta}
+\end{equation}
+Finally, we evaluate the order parameter term, using \eqref{eq:dFodeta} and \eqref{eq:d2Fodetadeta} which give
+\begin{equation}
+  \begin{aligned}
+    \frac{\delta^2F_\op[\eta_\star[\epsilon]]}{\delta\epsilon(x)\delta\epsilon(x')}
+    &=\frac{\delta}{\delta\epsilon(x)}\int dx''\,\frac{\delta\eta_\star[\epsilon](x'')}{\delta\epsilon(x')}\frac{\delta F_\op[\eta]}{\delta\eta(x'')}\bigg|_{\eta=\eta_\star[\epsilon]} \\
+    &=\int dx''\,\frac{\delta^2\eta_\star[\epsilon](x'')}{\delta\epsilon(x)\delta\epsilon(x')}\frac{\delta F_\op[\eta]}{\delta\eta(x'')}\bigg|_{\eta=\eta_\star[\epsilon]}
+    +\int dx''dx'''\frac{\delta\eta_\star[\epsilon](x'')}{\delta\epsilon(x)}\frac{\delta\eta_\star[\epsilon](x''')}{\delta\epsilon(x')}\frac{\delta^2F_\op[\eta]}{\delta\eta(x'')\delta\eta(x''')}\bigg|_{\eta=\eta_\star[\epsilon]} \\
+    &=\int dx''\,\frac{\delta^2\eta_\star[\epsilon](x'')}{\delta\epsilon(x)\delta\epsilon(x')}\big(b\epsilon(x)-g\epsilon(x)^2\eta_\star[\epsilon](x)\big)
+    +(b-2g\epsilon(x)\eta_\star[\epsilon](x))\frac{\delta\eta_\star[\epsilon](x)}{\delta\epsilon(x')} \\
+    &\qquad-g\int dx''\,\epsilon(x'')^2\frac{\delta\eta_\star[\epsilon](x'')}{\delta\epsilon(x)}\frac{\delta\eta_\star[\epsilon](x'')}{\delta\epsilon(x')}.
+  \end{aligned}
+\end{equation}
+Summing all three terms, we see a great deal of cancellation, with
+\[
+  \frac{\delta^2F[\eta_\star[\epsilon],\epsilon]}{\delta\epsilon(x)\delta\epsilon(x')}=C_0\delta(x-x')+g\eta_\star[\epsilon](x)^2\delta(x-x')-(b-2g\epsilon(x)\eta_\star[\epsilon](x))\frac{\delta\eta_\star[\epsilon](x)}{\delta\epsilon(x')}.
+\]
+We new need to evaluate this at $\langle\epsilon\rangle$. First, $\eta_\star[\langle\epsilon\rangle]=\langle\eta\rangle$, and
+\[
+  \frac{\delta^2F[\eta_\star[\epsilon],\epsilon]}{\delta\epsilon(x)\delta\epsilon(x')}\bigg|_{\epsilon=\langle\epsilon\rangle}=C_0\delta(x-x')+g\langle\eta(x)\rangle^2\delta(x-x')-(b-2g\langle\epsilon(x)\rangle\langle\eta(x)\rangle)\frac{\delta\eta_\star[\epsilon](x)}{\delta\epsilon(x')}\bigg|_{\epsilon=\langle\epsilon\rangle}.
+\]
+Computing the final functional derivative is the most challenging part. We will
+first compute its functional inverse, take the Fourier transform of that, and
+then use the basic relationship between Fourier functional inverses to find the
+form of the non-inverse. First, we note
+\begin{equation}
+  \frac{\delta^2F_\op[\eta]}{\delta\eta(x)\delta\eta(x')}\bigg|_{\eta=\langle\eta\rangle}
+  =\big[r-c_\perp\nabla_\perp^2-c_\parallel\nabla_\parallel^2+D\nabla_\perp^4+12u\langle\eta(x)\rangle^2\big]\delta(x-x'),
+\end{equation}
+which gives
+\begin{equation}
+  \begin{aligned}
+    \bigg(\frac{\delta\eta_\star[\epsilon](x)}{\delta\epsilon(x')}\bigg)^{\{-1\}}\bigg|_{\epsilon=\langle\epsilon\rangle}
+    &=\frac1{b-2g\langle\epsilon(x)\rangle\langle\eta(x)\rangle}\bigg[g\langle\epsilon(x)\rangle^2\delta(x-x')+\frac{\delta^2F_\op[\eta]}{\delta\eta(x)\delta\eta(x')}\bigg]_{\eta=\langle\eta\rangle} \\
+    &=\frac1{b-2g\langle\epsilon(x)\rangle\langle\eta(x)\rangle}\Big[
+      g\langle\epsilon(x)\rangle^2+r-c_\perp\nabla_\perp^2-c_\parallel\nabla_\parallel^2+D\nabla_\perp^4+12u\langle\eta(x)\rangle^2
+    \Big]\delta(x-x').
+  \end{aligned}
+\end{equation}
+Upon substitution of \eqref{eq:epsilon_star} and expansion to quadratic order it $\langle\eta(x)\rangle$, we find
+\begin{equation}
+  \begin{aligned}
+    \bigg(\frac{\delta\eta_\star[\epsilon](x)}{\delta\epsilon(x')}\bigg)^{\{-1\}}\bigg|_{\epsilon=\langle\epsilon\rangle}
+    &=\frac1b\Bigg\{r-c_\perp\nabla_\perp^2-c_\parallel\nabla_\parallel^2+D\nabla_\perp^4\\
+    &\qquad\qquad+\langle\eta(x)\rangle^2\bigg[12u+\frac{b^2g}{C_0^2}+\frac{2g}{C_0}(r-c_\perp\nabla_\perp^2-c_\parallel\nabla_\parallel^2+D\nabla_\perp^4)\bigg]+O(\langle\eta\rangle^4)\Bigg\}\delta(x-x').
+  \end{aligned}
+\end{equation}
+Defining $\widehat{\langle\eta\rangle^2}=\int dq'\,\langle\hat\eta(q')\rangle\langle\hat\eta(-q')\rangle$, its Fourier transform is then
+\begin{equation}
+  \begin{aligned}
+    G(q)
+    &=\frac1V\int dx\,dx'\,e^{-iq(x-x')}\bigg(\frac{\delta\eta_\star[\epsilon](x)}{\delta\epsilon(x')}\bigg)^{\{-1\}}\bigg|_{\epsilon=\langle\epsilon\rangle} \\
+    &=\frac1b\Bigg\{r+c_\perp q_\perp^2+c_\parallel q_\parallel^2+Dq_\perp^4+\widehat{\langle\eta\rangle^2}\bigg[12u+\frac{b^2g}{C_0^2}+\frac{2g}{C_0}(r+c_\perp q_\perp^2+c_\parallel q_\parallel^2+Dq_\perp^4)\bigg]+O(\langle\hat\eta\rangle^4)\Bigg\}.
+  \end{aligned}
+\end{equation}
+We can now compute $C(q)$ by taking its Fourier transform, using the convolution theorem for the second term:
+\begin{equation}
+  \begin{aligned}
+    C(q)
+    &=C_0+g\widehat{\langle\eta\rangle^2}-\int dq''\bigg(b\delta(q'')-\frac{gb}{C_0}\int dq'\langle\hat\eta_{q'}\rangle\langle\hat\eta_{q''-q'}\rangle\bigg)/G(q-q'') \\
+    &=C_0+g\widehat{\langle\eta\rangle^2}-b^2\bigg(\frac1{r+c_\perp q_\perp^2+c_\parallel q_\parallel^2+Dq_\perp^4}-\widehat{\langle\eta\rangle^2}\frac{12u+b^2g/C_0^2+\frac{2g}{C_0}(r+c_\perp q^2+c_\parallel q_\parallel^2+Dq_\perp^4)}{(r+c_\perp q_\perp^2+c_\parallel q_\parallel^2+Dq_\perp^4)^2}\bigg)\\
+    &\qquad+\frac{gb^2}{C_0}\int dq'\,dq''\frac{\langle\hat\eta_{q'}\rangle\langle\hat\eta_{q''-q'}\rangle}{r+c_\perp(q_\perp-q_\perp'')^2+c_\parallel(q_\parallel-q_\parallel'')^2+D(q_\perp-q_\perp'')^4}+O(\langle\hat\eta\rangle^4).
+  \end{aligned}
+\end{equation}
+Upon substitution of $\langle\hat\eta_q\rangle=\frac12\eta_*\big[\delta(q_\perp-q_*)+\delta(q_\perp+q_*)\big]\delta(q_\parallel)$, we have
+\begin{equation}
+  \begin{aligned}
+    C(q)
+    &=C_0+\frac14g\eta_*^2-b^2\bigg(\frac1{r+c_\perp q_\perp^2+c_\parallel q_\parallel^2+Dq_\perp^4}-\frac{\eta_*^2}4\frac{12u+b^2g/C_0^2+\frac{2g}{C_0}(r+c_\perp q^2+c_\parallel q_\parallel^2+Dq_\perp^4)}{(r+c_\perp q_\perp^2+c_\parallel q_\parallel^2+Dq_\perp^4)^2}\bigg)\\
+    &\qquad+\frac{gb^2\eta_*^2}{4C_0}\bigg(\frac2{r+c_\parallel q_\parallel^2+c_\perp q_\perp^2+Dq_\perp^4}+\frac1{r+c_\parallel q_\parallel^2+c_\perp(q_\perp-2q_*)^2+D(q_\perp-2q_*)^4} \\
+    &\qquad\qquad+\frac1{r+c_\parallel q_\parallel^2+c_\perp(q_\perp+2q_*)^2+D(q_\perp+2q_*)^4}\bigg)+O(\eta_*^4).
+  \end{aligned}
+\end{equation}
+Evaluating at $q=0$, we have
+\begin{equation}
+  \begin{aligned}
+    C(0)
+    &=C_0-\frac{b^2}r+\frac{\eta_*^2}4\bigg(g+\frac{b^2}{r^2}(12u+b^2g/C_0^2)+\frac{2gb^2}{C_0r}\frac{16Dq_*^4+3r}{8Dq_*^4+r}\bigg)
+  \end{aligned}
+  \label{eq:C0}
+\end{equation}
+Above the transition this has exactly the form of \eqref{eq:static_modulus} for
+any $g$; below the transition it has the same form at $g=0$ to order
+$\eta_*^2$. With $r=a\Delta T+c^2/4D+b^2/C_0$, $u=\tilde u-b^2g/2C_0^2$, and
+\begin{equation}
+  \eta_*^2=\begin{cases}
+    0 & \Delta T > 0 \\
+    -a\Delta T/3\tilde u & \Delta T \leq 0,
+  \end{cases}
+\end{equation}
+we can fit the ratios $b^2/a=1665\,\mathrm{GPa}\,\mathrm K$, $b^2/Dq_*^4=6.28\,\mathrm{GPa}$, and $b\sqrt{-g/\tilde u}=14.58\,\mathrm{GPa}$ with $C_0=(71.14-(0.010426\times T)/\mathrm K)\,\mathrm{GPa}$. The resulting fit the thin solid black line in Fig.~\ref{fig:data}.
+\end{widetext}
+
 \bibliographystyle{apsrev4-1}
 \bibliography{hidden_order}
 
-- 
cgit v1.2.3-70-g09d2