Random quench predicts universal properties of amorphous solids

At Kelvin-scale temperatures, glasses universally present mechanical and vibrational anomalies with respect to crystalline solids: below $1K$, the heat capacity behaves as $C(T)\propto T$, to be compared with $C(T)\propto T^{3}$ for crystals Anderson (1981). It is accepted that the anomalous behavior of $C(T)$ in glasses is caused by quantum mechanical tunneling between nearby energy minima in phase space, the two-level systems (TLS) initially proposed as a phenomenological model Anderson et al. (1972). Recently, a microscopic picture of the TLS has begun to emerge, thanks to numerical simulations in which quasi-localized vibrational modes have been identified and counted Lerner et al. (2016); Kapteijns et al. (2018); Lerner and Bouchbinder (2018); Angelani et al. (2018); Wang et al. (2019); Ji et al. (2019); Khomenko et al. (2020).

Moreover, near $10K$, glasses display an excess of vibrational modes over phonons, the so-called ’boson peak.’ It is not agreed what is the cause of the modes near the boson peak: the glass-specific behavior has variously been attributed to soft localized modes Buchenau et al. (1992); Parshin et al. (2007), generic stiffness disorder Schirmacher (2006); Schirmacher et al. (2007), proximity to jamming Wyart et al. (2005a, b); DeGiuli et al. (2014a), and proximity to elastic instability DeGiuli et al. (2014a).

The jamming approach predicts a regime in which the density of vibrational states $g(\omega)$ scales as $g(\omega)\propto\omega^{2}$ in any dimension $d$, sufficiently close to jamming and elastic instability DeGiuli et al. (2014a). The corresponding modes are extended but not plane waves, instead showing vortex-like motion at the particle scale. This law has been confirmed in numerical simulations DeGiuli et al. (2014a); Mizuno and Ikeda (2018); Shimada et al. (2018, 2020a), but is found to break down below some frequency $\omega_{0}$, below which there are only phonons and quasi-localized modes Mizuno and Ikeda (2018); Wang et al. (2019); Shimada et al. (2020a). The latter, now confirmed to be present in many glass models, has a density following $g(\omega)\propto\omega^{\alpha}$ where $3\leq\alpha\leq 4$ Lerner et al. (2016); Lerner and Bouchbinder (2017); Kapteijns et al. (2018); Lerner and Bouchbinder (2018); Angelani et al. (2018); Wang et al. (2019); Ji et al. (2019); Paoluzzi et al. (2019); Shimada et al. (2020a); Lerner (2020). Some authors argue that $\alpha=4$ in the thermodynamic limit Lerner (2020), while others argue that smaller exponents are possible due to interactions between localized instabilities Ji et al. (2019). Microscopic theory is needed to clarify these results.

The frequency $\omega_{0}$ setting the lower-limit of the jamming regime is controlled by the distance to elastic instability DeGiuli et al. (2014a). It was proposed that glasses dominated by repulsive interactions lay close to elastic instability due to the quench dynamics Wyart et al. (2005b); DeGiuli et al. (2014a). This suggests that a model faithful to the physics of the quench might shed light on the density of quasi-localized modes and the frequency $\omega_{0}$ below which they become important. Ideally, any such model should also reproduce the universal vibrational features captured in prevous models Schirmacher et al. (2007); DeGiuli et al. (2014a), such as the dip in sound speed and crossover in acoustic attenuation observed in many experiments Rufflé et al. (2006); Monaco and Giordano (2009); Baldi et al. (2011a); Ruta et al. (2012).

Here we present such a model. Following a crude but principled model of an overdamped quench into an inherent state (IS), we derive universal vibrational properties characterized by complex elastic moduli and the density of vibrational states. We recover the exact equations governing mean-field vibrational properties discussed previously DeGiuli et al. (2014a). As a bonus, our model predicts other universal mechanical features, namely short-range correlations in elastic moduli and long-range correlations in the IS stress, as observed recently in simulations Mizuno et al. (2016); Mizuno and Mossa (2019); Shakerpoor et al. (2020); Kapteijns et al. (2020); Lemaître (2015, 2017) and experiments ¹¹1These are inferred from strain measurements in colloidal glassesChikkadi et al. (2011); Jensen et al. (2014); Illing et al. (2016) and granular media Le Bouil et al. (2014).Wang et al. (2020). Ultimately, the model fails to predict the universal $\omega^{4}$ law discussed above, thus indicating the minimal features needed to recover mean-field results, while highlighting routes towards a more complete model.

Consider the elasticity equation

$\displaystyle\rho\frac{d^{2}\vec{u}_{i}}{dt^{2}}-\partial_{j}\left[S_{ijkl}\partial_{k}\vec{u}_{l}\right]=\vec{F}_{i},$

(1)

where $\rho$ is density, $\vec{u}_{i}$ is a displacement field, and $S_{ijkl}=C_{ijkl}+\delta_{ik}\sigma_{jl}$ in terms of the elastic modulus tensor $C_{ijkl}$ and the IS stress $\sigma_{jl}$. The elastic Green’s function $G_{ij}$ is the solution to (1) for a Dirac-delta function forcing, $\vec{F}_{j}(\vec{r})=\vec{f}_{j}\delta(\vec{r})$, that is, $\vec{u}_{i}(\vec{r})=G_{ij}(\vec{r})\cdot\vec{f}_{j}$. This is one of the fundamental observables in linear elasticity and captures many properties of linear response, including sound speed, damping due to the disorder, and the density of vibrational states, as discussed below.

At the mesoscale, the elastic moduli and the stress tensor $\sigma_{ij}$ can be considered to be spatially fluctuating fields. Vibrational properties can be derived from the disorder-averaged Green’s function $\overline{G}_{ij}$, for different models of random disorder. The model of Schirmacher corresponds to random Gaussian fluctuations of elastic moduli Schirmacher (2006); Schirmacher et al. (2007).

From the jamming approach, Ref. DeGiuli et al. (2014a) employed a microscopic lattice model, which does not directly correspond to (1). However, elastic instability is caused by the destabilizing effect of stress in particulate matter with short-range repulsive interactions Wyart et al. (2005b); DeGiuli et al. (2014a). Its proposed importance highlights stress as an important control parameter for vibrational properties. Moreover, it has been shown that the stress also plays a crucial role in the emergence of quasi-localized modes Lerner and Bouchbinder (2018). These works suggest that one should consider a model in which the IS stress $\sigma_{ij}$ is randomly fluctuating. Such an effort must immediately confront a no-go theorem of Di Donna and Lubensky DiDonna and Lubensky (2005). In a comprehensive treatment of non-affine correlations in random media, the latter authors showed that a random IS stress alone does not yield non-affine motion, and therefore cannot give rise to anomalous vibrational properties: a material that behaves affinely is a homogeneous continuum, the continuum limit of a crystal. How can we reconcile the importance of destabilizing stress with its apparently mild effect on non-affine motion?

We propose that the solution is to consider the quench itself. Indeed, as shown in DiDonna and Lubensky (2005), if random forces are added to an initially featureless continuum, then the relaxation to an IS will produce both a random IS stress and fluctuations in elastic moduli. The latter cause anomalous vibrational properties.

To probe how fluctuations in stress and elastic moduli are created during the final descent into an IS, we consider the following idealized model of a quench. We begin with the dense liquid in its natural state at the parent temperature $T_{0}$. This state is characterized by a stress tensor field $\tilde{\sigma}_{ij}(\vec{r})$, which we call the quench stress. In the case of a pair potential, this corresponds to (e.g. Morante et al. (2006))

$\displaystyle\tilde{\sigma}_{ij}(\vec{r})=\sum_{a}\delta(\vec{r}-\vec{r}_{a})\left[mv^{a}_{i}v^{a}_{j}+\mbox{$\frac{1}{2}$}\sum_{b\neq a}r^{ab}_{i}f^{ab}_{j}\right]$

(2)

where $\vec{v}^{a}$ is the velocity of particle $a$, $\vec{r}^{ab}=\vec{r}^{a}-\vec{r}^{b},$ and $\vec{f}^{ab}$ is the force exerted by particle $b$ on particle $a$.

From this state we then instantaneously quench the temperature to 0, so that the velocity field vanishes. Since the state is not in mechanical equilibrium, it will relax under the action of the quench stress $\tilde{\sigma}_{ij}$ until it finds a state of mechanical equilibrium, as depicted in Figure 1. The process ends as soon as a state of mechanical equilibrium is found. We will then compute the disorder-averaged Green’s function of the IS.

Since our aim is to model the emergence of structure in the glass, we wish to consider the simplest possible model of the dense liquid. Only the short-time response of the liquid is relevant, so we consider the latter as an initial homogeneous elastic continuum with elastic constants $\tilde{\lambda}$ and $\tilde{\mu}$, which gives the elastic modulus tensor $C_{ijkl}=\tilde{\lambda}\delta_{ij}\delta_{kl}+\tilde{\mu}(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})$. We assume, for simplicity, that these constants are spatially uniform. When the temperature is removed, the continuum deforms under the quench stress, just as a crumpled elastic membrane will relax when external constraints are removed. For simplicity we ignore inertia during the quench, so our model is one of an overdamped quench. This is equivalent to gradient descent as performed in numerical simulations.

Although crude, this model is perhaps the simplest that captures the final stages of descent into an IS; more advanced models could allow the liquid to have non-trivial structural correlations and consider inertia in the quenching process.

We aim to describe universal properties of this process. Since we work in the continuum, the relevant distribution of $\tilde{\sigma}_{ij}(\vec{r})$ can be inferred using field-theoretical arguments DeGiuli (2018a, b). Indeed, under our assumption of a featureless liquid, under coarse-graining the quench stress field will retain only short-range correlations. As we eventually seek the small $q$ behaviour of the solid, we take this small length scale to 0 immediately to obtain a simple Gaussian distribution of the symmetric tensor field $\tilde{\sigma}_{ij}(\vec{r})$:

$\displaystyle\mathbb{P}[\tilde{\sigma}]\propto\exp\left(-\mbox{$\frac{1}{2}$}\int_{r}\left[s_{1}\tilde{}\mathrlap{\!\not{\phantom{\sigma}}}\sigma_{ij}\tilde{}\mathrlap{\!\not{\phantom{\sigma}}}\sigma_{ij}+s_{2}\tilde{\sigma}_{ii}\tilde{\sigma}_{jj}\right]\right),$

(3)

where $\tilde{}\mathrlap{\!\not{\phantom{\sigma}}}\sigma_{ij}=\tilde{\sigma}_{ij}-\mbox{$\frac{1}{d}$}\delta_{ij}\tilde{\sigma}_{kk}$ is the deviatoric stress, and $s_{1}$ and $s_{2}$ are parameters related to the magnitude of quench stress by

	$\displaystyle\overline{\tilde{}\mathrlap{\!\not{\phantom{\sigma}}}\sigma_{ij}(\vec{r})\tilde{}\mathrlap{\!\not{\phantom{\sigma}}}\sigma_{ij}(0)}$	$\displaystyle=\frac{d^{2}-1}{s_{1}}\delta(\vec{r}),$		(4)
	$\displaystyle\overline{\tilde{\sigma}_{ii}(\vec{r})\tilde{\sigma}_{jj}(0)}$	$\displaystyle=\frac{1}{s_{2}}\delta(\vec{r})$		(5)

in $d$ dimensions. The constant component of the quench stress will be treated separately.

Standard renormalization arguments DeGiuli (2018a, b) predict that corrections to (3) will introduce a length scale $a$ which is on the order of the relevant microscopic length, namely the particle radius. We therefore expect our model to be valid on scales much larger than this size. In particular, the Gaussian distribution (3) should be valid, so long as the liquid state does not have relevant long-range correlations. In our continuum treatment, we will employ a cutoff $\Lambda\approx 2\pi/a$ in momentum space, so that corrections to (3) need not be explicitly incorporated.

Consider the displacement field $u_{i}(\vec{r})$ along the quench. At any moment, there is an instantaneous stress field $\sigma_{ij}[\vec{u}](\vec{r})$. In an overdamped quench, a new IS will be found as soon as $\sigma_{ij}$ is in mechanical equilibrium. (If inertia were considered, the system could jump over shallow minima before finding a resting state.) Di Donna and Lubensky found the new IS stress $\sigma_{ij}(\vec{r})$ and the elastic modulus tensor $C^{\prime}_{ijkl}(\vec{r})=C_{ijkl}+\delta C_{ijkl}(\vec{r})$ around the IS, to leading order in $u_{i}$, for arbitrary quench stress fields $\tilde{\sigma}_{ij}(\vec{r})$ with zero spatial average. The result is a pair of linear functionals

	$\displaystyle\delta C_{ijkl}(\vec{q})$	$\displaystyle=\mathscr{S}_{ijklmn}(\vec{q})\tilde{\sigma}_{mn}(\vec{q})$		(6)
	$\displaystyle\sigma_{ij}(\vec{q})$	$\displaystyle=\mathscr{P}_{ijkl}(\vec{q})\tilde{\sigma}_{kl}(\vec{q})$		(7)

in Fourier space, where $\mathscr{S}$ and $\mathscr{P}$ depend on the elastic moduli and the momentum $\vec{q}$. Explicitly, these take the form

	$\displaystyle\mathscr{S}_{ijklmn}(\vec{q})$	$\displaystyle=-\frac{1}{2(1-c)}[2c(\delta_{ij}\delta_{ke}\delta_{lf}+\delta_{kl}\delta_{ie}\delta_{jf})+(1-c)(\delta_{ik}\delta_{je}\delta_{lf}+\delta_{jl}\delta_{ie}\delta_{kf}+\delta_{il}\delta_{je}\delta_{kf}+\delta_{jk}\delta_{ie}\delta_{lf})]V_{efmn}$		(8)
	$\displaystyle\mathscr{P}_{ijkl}(\vec{q})$	$\displaystyle=\frac{1}{2}{P^{T}}_{ik}{P^{T}}_{jl}+\frac{1}{2}{P^{T}}_{il}{P^{T}}_{jk}-\frac{1}{1+2\mu}{P^{T}}_{ij}\hat{q}_{k}\hat{q}_{l},$		(9)

where $c=1/(1+2\mu)$, ${P^{T}}_{ij}=\delta_{ij}-\hat{q}_{i}\hat{q}_{j}$, and

$\displaystyle V_{efmn}=\hat{q}_{e}(\delta_{fm}\hat{q}_{n}+\delta_{fn}\hat{q}_{m})+\hat{q}_{f}(\delta_{em}\hat{q}_{n}+\delta_{en}\hat{q}_{m})-2(1+c)\hat{q}_{e}\hat{q}_{f}\hat{q}_{m}\hat{q}_{n}.$

(10)

Ref. DiDonna and Lubensky (2005) did not include any constant component of the quench stress, but this is easily added. To leading order in $u$, we find that if $\overline{\tilde{\sigma}_{ij}(\vec{r})}={\overline{p}}\delta_{ij}$, then this is equivalent to replacing the Lamé moduli by $\lambda=\tilde{\lambda}-d{\overline{p}}$ and $\mu=\tilde{\mu}+{\overline{p}}$.

We propose a model for universal properties of amorphous solids based on the quench into an inherent state. Under a single universal distribution of quench stress, our model predicts (i) short-range correlations of elastic moduli, as observed Mizuno et al. (2016); Kapteijns et al. (2020); Shakerpoor et al. (2020); (ii) long-range correlations of the IS stress, as observed Lemaître (2015, 2017); Mizuno and Mossa (2019); (iii) exact reduction to previous models, shown to reproduce universal vibrational anomalies DeGiuli et al. (2014a); Shimada et al. (2020b); and (iv) a tail of potentially unstable modes, beyond mean-field predictions, which leads to $g(\omega)\sim\omega^{3}$ in small systems and can rationalize larger exponents $g(\omega)\sim\omega^{\alpha}$ in large, stable systems.

New results include the relationship between the elastic-moduli fluctuations and the stress correlations; explicit expressions for the stress correlations in arbitrary dimension; justification of mean-field models which formerly were derived from lattice models; and the inevitable extension of mean-field models to have an eigenvalue tail in finite dimensions.

Our model is based on an overdamped quench, which enforces mechanical equilibrium in the IS, but does not enforce stability. This allows unstable configurations, like a pencil standing on its head. At the global level, stability is ensured by choosing parameters such that all vibrational modes have a positive real frequency; in particular we should have $e<e_{c}$ as described above. However, we do not enforce stability at each point in space. In particular, the Hessian field $H_{ij}(\vec{r},\vec{r}^{\prime})=\partial^{2}E/(\partial u_{i}(\vec{r})\partial u_{j}(\vec{r}^{\prime}))$ that controls local stability could have regions where it is not positive-definite. We predict Gaussian fluctuations of local elastic moduli, as observed in Mizuno et al. (2016); Shakerpoor et al. (2020), while it has very recently been argued that the moduli have a power-law tail due to localized modes Kapteijns et al. (2020). We expect that these modes are created by local relaxation in regions of an unstable Hessian. To rigorously predict the density of small-frequency localized modes in large systems, and their potential modifications to elastic moduli fluctuations, future work should thus add local stability as a remaining feature to the model.

In addition, we have assumed for simplicity that the initial pre-quenched liquid state has homogenous short-time elastic constants. It would be useful and relevant to allow these to have short-range Gaussian fluctuations, which will be transformed under the quench and may alter the elastic properties of the glass.

Acknowledgments: We are grateful to E. Lerner and M. Wyart for comments on the manuscript, and to the referees for useful comments. MS thanks H. Mizuno for useful information on simulations and A. Ikeda for his encouragement. MS is supported by JSPS KAKENHI Grant No. 19J20036.

Here we show how to derive the distribution of the inherent stress field. In particular we show that our model predicts the observed long-range stress correlations, and moreover connects the magnitude of these correlations to the coefficients controlling vibrational properties. We use $\overline{\cdot}$ to denote a disorder average, not to be confused with the Grassmann-conjugation defined elsewhere.

The distribution of $\sigma$ is

	$\displaystyle\mathbb{P}[\sigma]$	$\displaystyle=\overline{\prod_{ij}\delta\left[\sigma_{ij}-\mathscr{P}_{ijkl}\tilde{\sigma}_{kl}\right]}$
		$\displaystyle=\int{\mathscr{D}}\hat{\tau}\,\overline{\exp\left\{i\int_{q}\tau_{ij}(q)\left[\sigma_{ij}(q)-\mathscr{P}_{ijkl}(q)\tilde{\sigma}_{kl}(q)\right]\right\}}$
		$\displaystyle=\int{\mathscr{D}}\hat{\tau}\exp\left[i\int_{q}\tau_{ij}(q)\sigma_{ij}(q)\right]\overline{\exp\left[-i\int_{q}\tau_{ij}(q)\mathscr{P}_{ijkl}(q)\tilde{\sigma}_{kl}(q)\right]}$
		$\displaystyle=\int{\mathscr{D}}\hat{\tau}\exp\left[i\int_{q}\tau_{ij}(q)\sigma_{ij}(q)\right]\exp\left[-\int_{q}\tau_{ij}(q)\mathscr{P}^{\prime}_{ijkl}(q)\tau_{kl}(-q)\right],$

where

	$\displaystyle\mathscr{P}^{\prime}_{ijkl}(q)$	$\displaystyle=b_{1}\left(\frac{2\mu}{1+2\mu}\right)^{2}{P^{T}}_{ij}{P^{T}}_{kl}+\frac{b_{2}}{2}{P^{T}}_{ik}{P^{T}}_{jl}+\frac{b_{2}}{2}{P^{T}}_{il}{P^{T}}_{jk}+\frac{b_{2}}{(1+2\mu)^{2}}{P^{T}}_{ij}{P^{T}}_{kl}$
		$\displaystyle\equiv b_{1}^{\prime}{P^{T}}_{ij}{P^{T}}_{kl}+\frac{b_{2}}{2}\left({P^{T}}_{ik}{P^{T}}_{jl}+{P^{T}}_{il}{P^{T}}_{jk}\right).$

and the inverse tensor $K^{-1}_{klmn}=b_{1}\delta_{kl}\delta_{mn}+b_{2}\delta_{km}\delta_{ln}$ is determined by the relation $K_{ijkl}K^{-1}_{klmn}=\delta_{im}\delta_{jn}$. Then we have

	$\displaystyle b_{1}$	$\displaystyle=-\frac{1}{ds_{1}s_{2}}\left(s_{2}-\frac{s_{1}}{d}\right)$
	$\displaystyle b_{2}$	$\displaystyle=\frac{1}{s_{1}}.$

Since $\mathscr{P}^{\prime}_{ijkl}(q)q_{i}=\mathscr{P}^{\prime}_{ijkl}(q)q_{k}=0$, it is useful to perform a tensorial Helmholtz decomposition. We change the variable as $\tau_{ij}(q)=-iq_{i}\phi_{j}(q)-iq_{j}\phi_{i}(q)+\Psi_{ij}(q)$, where $q_{i}\Psi_{ij}(q)=q_{j}\Psi_{ij}(q)=0$. Thus, we have

	$\displaystyle\mathbb{P}[\sigma]$	$\displaystyle=\int{\mathscr{D}}\phi{\mathscr{D}}{\Psi}\exp\left\{i\int_{q}\left[-iq_{i}\phi_{j}(q)-iq_{j}\phi_{i}(q)+\Psi_{ij}(q)\right]\sigma_{ij}(q)-\int_{q}\Psi_{ij}(q)\mathscr{P}^{\prime}_{ijkl}(q)\Psi_{kl}(-q)\right\}$
		$\displaystyle=\int{\mathscr{D}}\phi{\mathscr{D}}{\Psi}\exp\left\{i\int_{r}\partial_{i}\phi_{j}(r)\left[\sigma_{ij}(r)+\sigma_{ji}(r)\right]-\int_{q}\left[\Psi_{ij}(q)\mathscr{P}^{\prime}_{ijkl}(q)\Psi_{kl}(-q)-i\Psi_{ij}(q)\sigma_{ij}(q)\right]\right\}.$

Decomposing ${\Psi}$ and $\sigma$ as ${\Psi}=({\Psi}+{\Psi}^{T})/2+({\Psi}-{\Psi}^{T})/2={\Psi}^{S}+{\Psi}^{A}$ and $\sigma=\sigma^{S}+\sigma^{A}$, we have

	$\displaystyle\mathbb{P}[\sigma]$	$\displaystyle=\delta\left[\sigma_{ij}(r)-\sigma_{ji}(r)\right]\delta\left[\partial_{i}\sigma_{ij}(r)\right]\int{\mathscr{D}}{\Psi}^{S}\exp\left\{-\int_{q}\left[\Psi_{ij}^{S}(q)\mathscr{P}^{\prime}_{ijkl}(q)\Psi_{kl}^{S}(-q)-i\Psi^{S}_{ij}(q)\sigma_{ij}(q)\right]\right\}$
		$\displaystyle\equiv\delta\left[\sigma_{ij}(r)-\sigma_{ji}(r)\right]\delta\left[\partial_{i}\sigma_{ij}(r)\right]\mathbb{P}^{S}[\sigma].$

The $\delta-$functions here are exactly the equations of mechanical equilibrium. The remaining symmetric part $\mathbb{P}^{S}[\sigma]$ further simplifies to

$\displaystyle\mathbb{P}^{S}[\sigma]$

$\displaystyle=\int{\mathscr{D}}{\Psi}^{S}\exp\left\{-\int_{q}\left[\Psi_{ij}^{S}(q)\left(b_{1}^{\prime}\delta_{ij}\delta_{kl}+b_{2}\delta_{ik}\delta_{jl}\right)\Psi_{kl}^{S}(-q)-i\Psi^{S}_{ij}(q)\sigma_{ij}(q)\right]\right\},$

where we used $q_{i}\Psi^{S}_{ij}(q)=0$.

In $d=2$, setting $\Psi^{S}_{ij}(q)=-\epsilon_{ik}\epsilon_{jl}q_{k}q_{l}\psi(q)$, we have

	$\displaystyle\mathbb{P}^{S}[\sigma]$	$\displaystyle=\int{\mathscr{D}}\psi\exp\left\{-\int_{q}\left[\epsilon_{im}\epsilon_{jn}\epsilon_{kp}\epsilon_{lq}\left(b_{1}^{\prime}\delta_{ij}\delta_{kl}+b_{2}\delta_{ik}\delta_{jl}\right)q_{m}q_{n}q_{p}q_{q}\psi(q)\psi(-q)+i\epsilon_{ik}\epsilon_{jl}q_{k}q_{l}\sigma_{ij}(q)\psi(q)\right]\right\}$
		$\displaystyle=\int{\mathscr{D}}\psi\exp\left\{-\int_{q}\left[(b_{1}^{\prime}+b_{2})q^{4}\psi(q)\psi(-q)+i\epsilon_{ik}\epsilon_{jl}q_{k}q_{l}\sigma_{ij}(q)\psi(q)\right]\right\}.$

To match the notation of DeGiuli (2018a, b), we set $1/2\tilde{\eta}=b_{1}^{\prime}+b_{2}$. Thus,

	$\displaystyle\mathbb{P}^{S}[\sigma]$	$\displaystyle=\int{\mathscr{D}}\psi\exp\left\{-\frac{\tilde{\eta}^{-1}}{2}\int_{q}q^{4}\left[\psi(q)+i\epsilon_{ik}\epsilon_{jl}\hat{q}_{k}\hat{q}_{l}q^{-2}\tilde{\eta}\sigma_{ij}(q)\right]\left[\psi(-q)+i\epsilon_{mp}\epsilon_{nq}\hat{q}_{p}\hat{q}_{q}q^{-2}\tilde{\eta}\sigma_{mn}(-q)\right]\right.$
		$\displaystyle\qquad\left.-\frac{\tilde{\eta}}{2}\int_{q}\epsilon_{ik}\epsilon_{jl}\hat{q}_{k}\hat{q}_{l}\sigma_{ij}(q)\epsilon_{mp}\epsilon_{nq}\hat{q}_{p}\hat{q}_{q}\sigma_{mn}(-q)\right\}$
		$\displaystyle=\frac{1}{Z_{\sigma}}\exp\left[-\frac{\tilde{\eta}}{2}\int_{q}\epsilon_{ik}\epsilon_{jl}\hat{q}_{k}\hat{q}_{l}\sigma_{ij}(q)\epsilon_{mp}\epsilon_{nq}\hat{q}_{p}\hat{q}_{q}\sigma_{mn}(-q)\right]$
		$\displaystyle=\frac{1}{Z_{\sigma}}\exp\left[-\frac{\tilde{\eta}}{2}\int_{r}\mbox{tr}^{2}\sigma(r)\right],$

where $Z_{\sigma}$ is the normalization constant and we used $q_{i}\sigma_{ij}(q)=0$. The coefficient $\tilde{\eta}$ is given by

$\displaystyle\tilde{\eta}$

$\displaystyle=\frac{(1+2\mu)^{2}}{2\{4\mu^{2}b_{1}+[1+(1+2\mu)^{2}]b_{2}\}}.$