# THz ultra-strong light-matter coupling up to 200K with continuously-graded
parabolic quantum wells

Paul Goulain  Centre de Nanosciences et de Nanotechnologies, CNRS UMR 9001,
University of Paris-Saclay, 91120 Palaiseau, France  Chris Deimert  Department
of Electrical and Computer Engineering, University of Waterloo, 200 University
Ave W, Waterloo, ON N2L 3G1, Canada  Mathieu Jeannin  Centre de Nanosciences
et de Nanotechnologies, CNRS UMR 9001, University of Paris-Saclay, 91120
Palaiseau, France  Stefano Pirotta  Centre de Nanosciences et de
Nanotechnologies, CNRS UMR 9001, University of Paris-Saclay, 91120 Palaiseau,
France  Wojciech Julian Pasek  Faculdade de Ciências Aplicadas, Universidade
Estadual de Campinas, Limeira, SP, 13484-350, Brazil  Zbigniew Wasilewski
Institute for Quantum Computing, University of Waterloo, 200 University Ave W,
Waterloo, ON N2L 3G1, Canada  Waterloo Institute for Nanotechnology,
University of Waterloo, 200 University Ave W, Waterloo, ON N2L 3G1, Canada
Department of Physics and Astronomy, University of Waterloo, 200 University
Ave W, Waterloo, ON N2L 3G1, Canada  Raffaele Colombelli
<EMAIL_ADDRESS>Centre de Nanosciences et de
Nanotechnologies, CNRS UMR 9001, University of Paris-Saclay, 91120 Palaiseau,
France  Jean-Michel Manceau<EMAIL_ADDRESS>Centre de Nanosciences et de Nanotechnologies, CNRS UMR 9001, University of
Paris-Saclay, 91120 Palaiseau, France

###### Abstract

Continuously graded parabolic quantum wells with excellent optical
performances are used to overcome the low-frequency and thermal limitations of
square quantum wells at terahertz frequencies. The formation of microcavity
intersubband polaritons at frequencies as low as 1.8 THz is demonstrated, with
a sustained ultra-strong coupling regime up to a temperature of 200 K. Thanks
to the excellent intersubband transition linewidth, polaritons present quality
factors up to 17. It is additionally shown that the ultra-strong coupling
regime is preserved when the active region is embedded in sub-wavelength
resonators, with an estimated relative strength
$\eta=\Omega_{R}/\omega_{0}=0.12$. This represents an important milestone for
future studies of quantum vacuum radiation because such resonators can be
optically modulated at ultrafast rates, possibly leading to the generation of
non-classical light via the dynamic Casimir effect. Finally, with an effective
volume of $2.10^{-6}\lambda_{0}^{3}$, it is estimated that fewer than 3000
electrons per resonator are ultra-strongly coupled to the quantized
electromagnetic mode, proving it is also a promising approach to explore few-
electron polaritonic systems operating at relatively high temperatures.

Polariton, Intersubband, Metamaterial, Terahertz

##  I Introduction

Following the pioneering work by Haroche et al.[1], engineering of the light-
matter interaction has been one of the most active fields of modern physics.
Remarkably, when the losses in the system are lower than the coupling
strength, a new regime is attained - the strong coupling regime (SCR) - where
light and matter exchange energy coherently and periodically. In the frequency
domain, this leads to a radical change of the system’s spectral response. From
the first observation with Rydberg atoms, the SCR has been demonstrated in a
plethora of systems spanning excitons, organic molecules, electronic
transitions, superconducting qubits and many others.[2] The strength of the
light-matter interaction is often gauged by a dimensionless parameter $\eta$
that is the ratio between the coupling constant $\Omega_{R}$ (also called the
vacuum Rabi frequency) over the resonant transition frequency $\omega_{0}$.[3]
Above a value $\mathrm{\eta>0.1}$, one enters the ultra-strong coupling (USC)
regime where the diamagnetic terms of the interaction Hamiltonian start to
play an important role, leading to a deviation from the linear approximation
and the formation of a sizeable population of virtual photons in the ground
state of the system.[4] The same foundational article by Ciuti et al.[4] also
proposed that an abrupt modulation of the system ground state leads to a
release of such virtual population as real photons, an approach that could
lead to the development of non-classical light emitters at long-wavelengths.

The first candidate for this USC regime of interaction were intersubband (ISB)
polaritons.[4] These hybrid states stem from the strong coupling between an
ISB transition (more rigorously, an ISB plasmon[5]) within the conduction band
of doped semiconductor quantum wells and the quantized electro-magnetic mode
of a microcavity. They can operate down to the THz frequency range, while
their Rabi frequency - proportional to the square root of the introduced
dopant density - can reach values equivalent to the bare transition. Soon
after the first experimental demonstration[6], the USC regime with ISB
polaritons was demonstrated in two different microcavity configurations.[7, 8]
Since then, USC was demonstrated with organic systems operating at visible
wavelength and ambient temperature as thoroughly discussed in [3], and more
recently with plasmonic resonances embedded in microcavity that could be in
principle extended to the far infrared range.[9] With doped heterostructures,
the coupling strength was increased further using two different strategies
that alter the character of the quantum well (QW) resonance. In the first
configuration, a massive increase of the doping density in QWs was found to
lead to the formation of a multisubband plasmon[10]: a resonance essentially
locked to the plasma frequency. This enabled a coupling strength
$\mathrm{\eta=0.47}$, once placed in a resonant metallic microcavity at mid-
infrared wavelengths.[11] In the second configuration, the application of a
strong magnetic field to a 2D electron gas leads to the formation of Landau
levels at sub-THz frequencies. When coupled to sub-wavelength resonators,
Landau polaritons form with a coupling strength that can reach values beyond
unity.[12, 13, 14] Several attempts have been made to obtain ultrafast
modulation of the USC regime with ISB polaritons,[15, 16] with most of the
activity aiming at measuring quantum vacuum radiation now conducted with
Landau polaritons.[17] Although they offer coupling strengths around/beyond
unity, they require elevated magnetic fields and cryogenic temperatures below
8 K.

Two primary factors have prevented further increase of $\eta$ with THz ISB
polaritons - particularly in square QWs. On one hand, any attempt to increase
the coupling strength ($\Omega_{R}$) by increasing the doping results in a
blueshift of the ISB transition, leaving the ratio $\eta$ unchanged. This is
due to the Coulomb interactions that act as a self-energy correction on the
optical transition. This effect, commonly known as the depolarization
shift,[5, 18] is proportional to the square root of the sheet doping density.
On the other hand, any attempt to lower the transition frequency
($\omega_{0}$) while maintaining a high doping density will spread the
population ($\mathrm{n_{2D}}$) over two subbands. This leads to the activation
of a second ISB transition and no further optimization of $\eta$ is possible.
Careful simulations of a square QW with a doping density of 1.1011 cm-2, using
a commercial solver (NextNano), reveal that it is not possible to obtain a
single and sharp ISB transition below 2.5 THz. Accessing the lower frequency
range is, however, a prerequisite to demonstrate non-classical light
generation via the dynamic Casimir effect. The reason is that elucidating such
subtle effects requires coherent field correlation measurements, which have so
far been done using the electro-optic sampling technique,[19, 20, 21, 22]
whose highest sensitivity bandwidth is in the 1-2 THz range.

##  II Quantum wells with parabolic energy potential

###  II.1 Design and implementation

A promising approach to circumvent this problem is to move from the standard
square QW to parabolic energy potentials. The key advantage is that the
electron-electron interaction cancels out - analogous to the Kohn’s theorem
for cyclotron resonance[23] \- leaving the transition energy independent from
distribution and density of the electrons. Digitally graded parabolic energy
potentials have already been used for the demonstration of strong light-matter
coupling at THz frequencies, but they suffer low optical performances due to
the interface roughness broadening that is inherent to the interdigitated
epitaxial growth technique.[24, 25] Furthermore, these demonstrations were not
designed to operate below the frequency of 2.5 THz imposed by classical square
QW structures. A recent experimental tour de force has radically changed the
perspective in that field. In Ref.[26], the parabolic energy potential is no
longer mimicked by a series of narrow square QWs, but it is faithfully
reproduced by a continuously graded alloy composition. This requires precise
time-varying group III fluxes, which are challenging to produce in a
conventional MBE setup, primarily because there is a thermal lag in the
response of the group III cell. The approach developed by Deimert at al. in
Ref.[26] uses a transfer function model to counteract the thermal dynamics,
along with an additional corrective step to achieve high-precision composition
gradients at typical growth rates. This approach led to high performance THz
ISB transitions with record linewidths up to 150 K, and to stable operation up
to room temperature.[27] In the present letter, we build on that result and we
experimentally show USC regime with ISB transitions below the inherent
frequency limit of square QWs.

Figure 1:  (a) Schrödinger-Poisson simulation of the parabolic quantum well
with a continuously graded alloy. In black is the Conduction band energy, the
green lines are the squared moduli of the wave functions
$\left|\Psi\right|^{2}$, while the dash line marks the energy position of the
Fermi level. (b) Calculated absorption of the system with a linewidth of 10%.

The heterostructure comprises a stack of 4 PQWs with 1.1011 cm-2 nominal sheet
electron density per well. The wells were realized in the AlxGa1-xAs material
system by varying the aluminum composition x in the range 2%-20% along the
growth direction. The PQWs were each 130.75 nm wide, separated by 20 nm
Al0.2Ga0.8As barriers where the Si doping is introduced as $\delta$-layers at
the barrier center. Two GaAs cap layers (50-nm thick) are placed on top and
bottom of the structure. The overall active region thickness is 783 nm. Figure
1a shows the Schrödinger-Poisson simulation of the grown structure performed
with the commercial software nextnano[28] and the calculated absorption
centered at 2.2 THz (Figure 1b).

###  II.2 Characterization of the optical performances

We have first characterized the optical performance of the grown structure.
The sample is shaped in a multi-pass prism configuration with the facets
polished at 45 degrees angle. Ti/Au surface coating is made on the side of the
epitaxial growth to force the boundary conditions at the semiconductor-metal
interface and ensure a maximum of electric field on the parabolic quantum well
repetition. The waveguides were placed into a continuous-flow cryostat inside
a Fourier transform infrared spectrometer (Bruker IFS66V), with polarized
broadband THz light incident on the facet. To reach a spectral resolution as
low as 0.5 cm-1, the whole FTIR was placed under vacuum to remove the
contribution of water absorption lines that are present even at low relative
humidity levels $\mathrm{(RH<1\%)}$. Detection was performed with a liquid-He-
cooled Si bolometer. Figure 2a shows the color-map evolution of the
transmittance as function of the temperature. By modeling the ISB transition
with a diagonal dielectric tensor as described in Refs.[29, 30, 31], the
central frequency ($\nu_{0}$), the linewidth ($\gamma_{i}$) and the plasma
frequency ($\omega_{p}$) can be extracted. Furthermore, one can gauge the
evolution of the stack absorption as a function of the temperature. For a
single well, the integrated 2D absorption is expressed as in [32]:

$\int_{-\infty}^{+\infty}\alpha_{2D}(\omega)d\omega=\frac{\omega^{2}_{p}L_{QW}}{c\sqrt{\varepsilon_{s}}}$
(1)

Hence, one can define the parameter $\alpha=\sqrt{N_{QW}}\omega_{p}$, where
$N_{QW}$ is the number of wells, and monitor the stack absorption strength as
a function of the temperature.

The evolution of the linewidth ($\gamma_{i}$) and the central frequency
($\nu_{0}$) are plotted in Figure 2b and c respectively. We first note a very
stable single absorption peak at a frequency around 2.1 THz over the whole
temperature range. The linewidth shows stronger temperature dependence, with a
record value of 69 GHz at 8 K (3.2% FWHM) that remains below 10% up to 150 K.
Despite having several subbands populated at higher temperature, the parameter
alpha (Figure 2d) is stable over the whole temperature range, evidence that
the number of electrons participating in the absorption remains mostly
constant thanks to the harmonicity of the system.

Figure 2:  (a) Transmittance measurement of the grown heterostructure as
function of the temperature. (b) Extracted linewidth using Voigt fitting
formula as function of the temperature. We measured 69 GHz linewidth at 10 K.
(c) Central frequency of the transition as function of the temperature. The
dashed line is a guide for the eye to show the stability of the transition
central frequency (d) Integrated absorption of the transition as function of
the temperature.

##  III Ultra-Strong light-matter coupling in microcavity

Figure 3:  (a) Experimental reflectance of the ISB polaritonic system as
function of the temperature. The dots mark the positions of the reflectance
minima (b) Fit (dashed line) of the experimental reflectance at 78 K (full
line), using finite element simulations. In inset is a schematic of the micro-
ribbon resonator.(c) Simulation of the polaritonic system dispersion,
attesting that the system is at the anti-crossing. The dots mark the
experimental positions of the polaritonic peaks in the sample having a ribbon
size of 18 µm. (d) Coupling strength parameter ($\eta=\Omega_{R}/\omega_{0}$)
as function of the temperature. The dashed line marks the limit of the ultra-
strong coupling regime.

We then explore operation in the strong light-matter coupling regime. To this
scope, the active region is placed within non-dispersive metal-insulator-metal
(MIM) ribbon micro-resonators.[33] The details of the micro-fabrication
process can be found elsewhere.[34] We have fabricated several samples with
different sizes of metallic ribbon that we have measured in reflectance at 300
K. For each samples, the minima of reflectance from the different resonant
photonic modes ($\mathrm{TM_{01},TM_{03}}$) were placed on the dispersion plot
obtained from a full EM simulation along with the experimentally measured
position of the ISB transition. Hence, we could define the sample placed at
the anti-crossing that allows measuring the Rabi splitting and accurately
estimate the Rabi frequency. The lateral size of the metal stripe, $\Lambda$,
is 18 µm, which sets the resonator operating frequency. The gap in between
each stripe, D, is equal to 4 µm. Figure 3a shows the color-map evolution of
the reflectance as function of the temperature. One can clearly see the lift
of degeneracy between cavity and ISB transition that occurs already below 250
K, and the two polaritonic branches gaining in contrast while the temperature
decreases. To infer the Rabi frequency of the strongly coupled system, we have
performed careful EM simulations to fit the measured Rabi splitting, where the
ISB transition is depicted using the diagonal dielectric tensor approach
aforementioned. In the quantum well plane, the material is described by the
GaAs permittivity with phonons and a Drude model accounting for the free
motion of the electrons. Doing so, the only fitting parameter is the sheet
doping density ($\mathrm{n_{2D}}$) that enters the plasma frequency
expression:

$\omega_{p}=\sqrt{\frac{f_{12}e^{2}n_{2D}}{\varepsilon_{0}\varepsilon_{r}m^{*}L_{QW}}}$
(2)

Here $f_{12}$ is the oscillator strength (equal to 0.97 for a PQW) between the
first and second subband, $m^{*}$ is the GaAs effective mass and $L_{QW}$ is
the length of the period comprising one well and one barrier. Figure 3b shows
the experimental fitted reflectance at 10 K, where the $n_{2D}$ is equal to
4.1010 cm-2, a factor two below the nominal doping. This discrepancy with the
nominal doping value might be due to the depletion of the two QWs close to the
metal-semiconductor interface where strong band bending occurs. Based on the
fitted permittivity, we ran a full simulation with different metal stripe
sizes in order to reconstruct the dispersion and confirm that our microcavity
is indeed at the anti-crossing condition (Figure 3c). It is important to note
that experimentally, we have access to the Rabi splitting frequency, expressed
as $\Omega_{\mathrm{splitting}}=2\Omega_{R}\sqrt{\Gamma_{\mathrm{opt}}}$.
Here, $\mathrm{\Gamma_{opt}}$ is a dimensionless parameter, which represents
the fraction of electromagnetic energy coupled to the out-of-growth-plane
component of the electric field and spatially overlapping the active region,
as described in Ref.[35]. In turn, isolating the Rabi frequency in the above
expression, and using $\mathrm{\Gamma_{opt}}$=0.95 (obtained from the finite
element method simulations), we estimate that $\Omega_{R}$ = 0.25 THz. We can
extract $\Omega_{R}$ along with $\omega_{0}$ for the whole set of temperatures
and estimate the dimensionless parameter that gauges the coupling strength
$\eta=\Omega_{R}/\omega_{0}$. Figure 3d demonstrates that the system operates
in the ultra-strong coupling regime up to 200 K. Furthermore, we have
extracted the linewidth of the polaritonic states, showing $\gamma_{pol}$ of
around 180 GHz and 140 GHz up to 78K (for lower and upper polaritonic states
respectively). This is equivalent to Q factors of 10 and 17, which is about 5
times higher than the one of Landau polaritons. Such performance is
particularly important in view of non-adiabatic switching experiments since
the emission rate of vacuum quanta can be significantly increased by a
narrower polariton linewidth, as thoroughly discussed in. [36] We also
anticipate that the increase in the transition coherence time shall improve
the signal to noise ratio within the coherent electro-optic sampling
technique, in comparison to the incoherent background radiation inherently
present.[19, 20, 21, 22] Another way to gauge the importance of improved
polaritons linewidth is to use the figure of merit defined in [37], that is
$U=\sqrt{\eta C}$, where $\eta$ is the coupling strength and C the
cooperativity defined as $C=4\Omega_{R}^{2}/\kappa\gamma_{i}$, where $\kappa$
and $\gamma_{i}$ are the cavity and transition linewidth respectively. A value
of U$\geq$1 indicates that the system operates in the USC regime with a high
degree of coherence, allowing to access its physics without being hindered by
decoherence. For the present work, even though the coupling strength is
relatively small, we estimate U equals to 1.7, thanks to the improved
transition linewidth. In our future work, we will aim at improving the
coupling strength by lowering the transition to 1.1 THz while increasing the
doping density by a reasonable factor 3 ($\sqrt{3}$ on $\Omega_{R}$). This
would lead to an increase of $\eta$ by a factor 3.5 and in turn to U = 6.7 if
the ISB transition linewidth is preserved. Such value would in fact be
equivalent to the one reported in the literature for Landau polaritons,
attesting that our material system is an alternative to it.

##  IV Ultra-Strong light matter coupling in sub-wavelength resonators

Figure 4:  (a) Scanning electron microscope image of the LC resonators array.
Colorized in red is the capacitive section hosting the parabolic quantum
wells, the rest being gold (b) Simulated confined electric field in the
direction orthogonal to the growth plane (Ez). (c) Dispersion of the
polaritonic system as a function of the antenna length. Dashed lines are the
secular equation while the dots are the experimental points. The dotted
parallel lines mark the edges of the polaritonic gap. (d) Experimental
reflectance of the array with antenna length of 11 µm at a temperature of 78
K, attesting of the excellent optical performances of the system.

Finally, we have verified the compatibility of this active region with
recently developed deeply sub-wavelength microcavities. The idea is to use a
cavity that operates in the high-frequency inductor-capacitor (LC) resonant
mode and which confines the EM field within a sub-wavelength volume.[38, 39]
Such resonators have been successfully used in recent years to demonstrate
strong-light matter coupling with a reduced number of electrons [40] and to
increase the operating temperature of QWIPs.[41] More recently, in pursuit of
non-adiabatic cavity modulation, the implementation of a frequency switching
functionality of 280 GHz in less than 200 fs was demonstrated by taking
advantage of the circuital properties of the device.[42] Figure 4a shows a
scanning electron microscope (SEM) picture of the fabricated device. The
fabrication details can be found in Ref.[38]. The capacitor section (colorized
red) is a cylinder of 1.5 µm diameter by 783 nm height, with an effective
volume of $2.10^{-6}\lambda_{0}^{3}$. Since the resonance frequency scales
with $\sqrt{LC}$, we have fabricated several samples with different antenna
lengths to reconstruct the polaritonic dispersion (Figure 4b). From the
reflectivity minima, we can identify the two polariton branches and fit them
with the secular equation[8]

$\left(\omega-\omega_{c}^{2}\right)\left(\omega-\tilde{\omega}_{21}^{2}\right)=\Gamma_{\mathrm{opt}}\omega_{p}^{2}\omega_{c}^{2}$
(3)

Where $\omega_{c}$ is the cavity resonance numerically calculated,
$\tilde{\omega}_{21}^{2}$ is the measured transition while
$\Gamma_{\mathrm{opt}}$ and $\omega_{\mathrm{p}}$ were defined above. The
dotted lines are placed along the asymptotes and reveal the presence of a
polaritonic gap, a key signature of the USC regime. Figure 4c is the
reflectance measurement at the anti-crossing point where we measure a Rabi
splitting $\Omega_{s}$ = 0.37 THz. The optical confinement within the LC
resonator is lower than the MIM cavities, $\Gamma_{\mathrm{opt}}$ = 0.83 in
this case. We estimate a Rabi frequency of 0.24 THz and a coupling strength
$\eta$ of 0.11, confirming that the system operates in the USC regime. Since
the electron sheet density is $\mathrm{n_{2D}=4.10^{10}}$ cm-2, we estimate
that the USC is achieved with around 3000 electrons per resonator coupled to
the quantized EM mode. Even lower numbers could be reached by reducing the
size of the device capacitance. For instance, considering a single PQW active
region and leveraging the higher resolution offered by the electron beam
lithography, a resonator with a capacitive section of 500 nm diameter over 315
nm height could be fabricated. With such a device, we estimate that the USC
could be achieved with only 400 electrons per resonator, a further step toward
the study of the USC with few electrons.[43, 44]

##  V Conclusion

In conclusion, we have used continuously graded parabolic QWs with outstanding
optical performances to demonstrate ISB polaritons around 2 THz. The intrinsic
low frequency bottleneck of SQWs is overcome thanks to the robustness of the
harmonic oscillator allowing the formation of a polariton mode as low as 1.8
THz with Q factors up to 17. The system operates in the ultra-strong light-
matter coupling regime up to 200 K. We further demonstrate that the USC regime
is maintained in deeply sub-wavelength volumes where around 3000 electrons per
resonator are involved in the interaction. A further reduction towards
hundreds of electrons could be envisioned using state-of-the-art
nanofabrication processing. Finally, further optimization of the system
cooperativity can be envisioned along with the implementation of an ultrafast
switching element as in Ref.[42] that would make such system a more practical
alternative to Landau polaritons in the study of quantum vacuum radiation at
long-wavelengths.

Acknowledgements

This work was supported by the European Union Future and Emerging Technologies
(FET) Grant No. 737017 (MIR-BOSE), and from the French National Research
Agency, project TERASEL (ANR-18-CE24-0013). This work was partially supported
by the French RENATECH network, Canada First Research Excellence Fund, and the
Natural Sciences and Engineering Research Council of Canada (NSERC). We thank
C. Ciuti for granting access to Paris Cité University cleanroom, and P.
Filloux and S. Suffit for help.

## References

  * [1] Y. Kaluzny, P. Goy, M. Gross, J. M. Raimond, S. Haroche,  _Physical Review Letters_ 1983, _51_ , 13 1175.
  * [2] D. N. Basov, A. Asenjo-Garcia, P. J. Schuck, X. Zhu, A. Rubio,  _Nanophotonics_ 2021, _10_ , 1 549.
  * [3] A. Frisk Kockum, A. Miranowicz, S. De Liberato, S. Savasta, F. Nori,  _Nature Reviews Physics_ 2019, _1_ , 1 19.
  * [4] C. Ciuti, G. Bastard, I. Carusotto,  _Physical Review B_ 2005, _72_ , 11.
  * [5] T. Ando, A. B. Fowler, F. Stern,  _Reviews of Modern Physics_ 1982, _54_ , 2 437.
  * [6] D. Dini, R. Kohler, A. Tredicucci, G. Biasiol, L. Sorba,  _Physical Review Letters_ 2003, _90_ , 11.
  * [7] A. A. Anappara, S. De Liberato, A. Tredicucci, C. Ciuti, G. Biasiol, L. Sorba, F. Beltram,  _Physical Review B_ 2009, _79_ , 20 201303.
  * [8] Y. Todorov, A. M. Andrews, R. Colombelli, S. De Liberato, C. Ciuti, P. Klang, G. Strasser, C. Sirtori,  _Physical Review Letters_ 2010, _105_ , 19 196402.
  * [9] D. G. Baranov, B. Munkhbat, E. Zhukova, A. Bisht, A. Canales, B. Rousseaux, G. Johansson, T. J. Antosiewicz, T. Shegai,  _Nature Communications_ 2020, _11_ , 1 2715.
  * [10] F. Alpeggiani, L. C. Andreani,  _Physical Review B_ 2014, _90_ , 11.
  * [11] B. Askenazi, A. Vasanelli, Y. Todorov, E. Sakat, J.-J. Greffet, G. Beaudoin, I. Sagnes, C. Sirtori,  _ACS Photonics_ 2017, _4_ , 10 2550.
  * [12] V. M. Muravev, I. V. Andreev, I. V. Kukushkin, S. Schmult, W. Dietsche,  _Physical Review B_ 2011, _83_ , 7 075309.
  * [13] G. Scalari, Maissen C., Turčinková D., Hagenmüller D., De Liberato S., Ciuti C., Reichl C., Schuh D., Wegscheider W., Beck M., Faist J.,  _Science (New York, N.Y.)_ 2012, _335_ , 6074 1323.
  * [14] A. Bayer, M. Pozimski, S. Schambeck, D. Schuh, R. Huber, D. Bougeard, C. Lange,  _Nano Letters_ 2017, _17_ , 10 6340.
  * [15] G. Günter, A. A. Anappara, J. Hees, A. Sell, G. Biasiol, L. Sorba, S. De Liberato, C. Ciuti, A. Tredicucci, A. Leitenstorfer, R. Huber,  _Nature_ 2009, _458_ , 7235 178.
  * [16] S. Zanotto, R. Degl’Innocenti, J.-H. Xu, L. Sorba, A. Tredicucci, G. Biasiol,  _Physical Review B_ 2012, _86_ , 20 201302.
  * [17] M. Halbhuber, J. Mornhinweg, V. Zeller, C. Ciuti, D. Bougeard, R. Huber, C. Lange,  _Nature Photonics_ 2020, _14_ , 11 675.
  * [18] M. Załużny,  _Physica Status Solidi_ 1982, _113_ , 1 K1.
  * [19] A. S. Moskalenko, C. Riek, D. V. Seletskiy, G. Burkard, A. Leitenstorfer,  _Physical Review Letters_ 2015, _115_ , 26 263601.
  * [20] C. Riek, D. V. Seletskiy, A. S. Moskalenko, J. F. Schmidt, P. Krauspe, S. Eckart, S. Eggert, G. Burkard, A. Leitenstorfer,  _Science (New York, N.Y.)_ 2015, _350_ , 6259 420.
  * [21] I.-C. Benea-Chelmus, C. Bonzon, C. Maissen, G. Scalari, M. Beck, J. Faist,  _Physical Review A: Atomic, Molecular, and Optical Physics_ 2016, _93_ , 4 043812.
  * [22] I.-C. Benea-Chelmus, F. F. Settembrini, G. Scalari, J. Faist,  _Nature_ 2019, _568_ , 7751 202.
  * [23] W. Kohn,  _Physical Review_ 1961, _123_ , 4 1242.
  * [24] M. Geiser, F. Castellano, G. Scalari, M. Beck, L. Nevou, J. Faist,  _Physical Review Letters_ 2012, _108_ , 10 106402.
  * [25] B. Paulillo, J.-M. Manceau, L. H. Li, A. G. Davies, E. H. Linfield, R. Colombelli,  _Applied Physics Letters_ 2016, _108_ , 10 101101.
  * [26] C. Deimert, Z. R. Wasilewski,  _Journal of Crystal Growth_ 2019, _514_ 103\.
  * [27] C. Deimert, P. Goulain, J.-M. Manceau, W. Pasek, T. Yoon, A. Bousseksou, N. Y. Kim, R. Colombelli, Z. R. Wasilewski,  _Physical Review Letters_ 2020, _125_ , 9.
  * [28] S. Birner, T. Zibold, T. Andlauer, T. Kubis, M. Sabathil, A. Trellakis, P. Vogl,  _IEEE Transactions on Electron Devices_ 2007, _54_ , 9 2137.
  * [29] M. Załużny, C. Nalewajko,  _Physical Review B_ 1999, _59_ , 20 13043.
  * [30] J. M. Manceau, S. Zanotto, T. Ongarello, L. Sorba, A. Tredicucci, G. Biasiol, R. Colombelli,  _Applied Physics Letters_ 2014, _105_ , 8.
  * [31] M. Knorr, J. M. Manceau, J. Mornhinweg, J. Nespolo, G. Biasiol, N. L. Tran, M. Malerba, P. Goulain, X. Lafosse, M. Jeannin, M. Stefinger, I. Carusotto, C. Lange, R. Colombelli, R. Huber,  _Physical Review Letters_ 2022, _128_ , 24 247401.
  * [32] M. Helm,  In H. C. Liu, F. Capasso, editors, _Intersubband Transitions in Quantum Wells: Physics and Device Applications I_ , volume 62 of _Semiconductors and Semimetals_ , 1–99. Elsevier, 1999.
  * [33] Y. Todorov, L. Tosetto, J. Teissier, A. M. Andrews, P. Klang, R. Colombelli, I. Sagnes, G. Strasser, C. Sirtori,  _Optics Express_ 2010, _18_ , 13 13886.
  * [34] J. M. Manceau, S. Zanotto, I. Sagnes, G. Beaudoin, R. Colombelli,  _Applied Physics Letters_ 2013, _103_ , 9.
  * [35] S. Zanotto, R. Degl’Innocenti, L. Sorba, A. Tredicucci, G. Biasiol,  _Physical Review B_ 2012, _85_ , 3 035307.
  * [36] S. De Liberato, C. Ciuti, I. Carusotto,  _Physical Review Letters_ 2007, _98_ , 10.
  * [37] P. Forn-Díaz, L. Lamata, E. Rico, J. Kono, E. Solano,  _Reviews of Modern Physics_ 2019, _91_ , 2 025005.
  * [38] B. Paulillo, J. M. Manceau, A. Degiron, N. Zerounian, G. Beaudoin, I. Sagnes, R. Colombelli,  _Optics Express_ 2014, _22_ , 18 21302.
  * [39] C. Walther, G. Scalari, M. I. Amanti, M. Beck, J. Faist,  _Science (New York, N.Y.)_ 2010, _327_ , 5972 1495.
  * [40] M. Jeannin, G. Mariotti Nesurini, S. Suffit, D. Gacemi, A. Vasanelli, L. Li, A. G. Davies, E. Linfield, C. Sirtori, Y. Todorov,  _ACS Photonics_ 2019, _6_ , 5 1207.
  * [41] M. Jeannin, T. Bonazzi, D. Gacemi, A. Vasanelli, S. Suffit, L. Li, A. G. Davies, E. Linfield, C. Sirtori, Y. Todorov,  _Applied Physics Letters_ 2020, _117_ , 25 251102.
  * [42] P. Goulain, A. D. Koulouklidis, J.-M. Manceau, C. Daskalaki, B. Paulillo, K. Maussang, S. Dhillon, J. R. Freeman, L. Li, E. H. Linfield, S. Tzortzakis, R. Colombelli,  _ACS Photonics_ 2021, _8_ , 4 1097.
  * [43] Y. Todorov, C. Sirtori,  _Physical Review X_ 2014, _4_ , 4 041031.
  * [44] J. Keller, G. Scalari, S. Cibella, C. Maissen, F. Appugliese, E. Giovine, R. Leoni, M. Beck, J. Faist,  _Nano Letters_ 2017, _17_ , 12 7410.


# Renormalization of shell model of turbulence

Mahendra K. Verma<EMAIL_ADDRESS>Department of Physics, Indian Institute of
Technology Kanpur, Kanpur 208016, India  Shadab Alam<EMAIL_ADDRESS>Department of Mechanical Engineering, Indian Institute of Technology Kanpur,
Kanpur 208016, India

###### Abstract

Renormalization enables a systematic scale-by-scale analysis of multiscale
systems. In this paper, we employ renormalization group (RG) to the shell
model of turbulence and show that the RG equation is satisfied by
$|u_{n}|^{2}=K_{\mathrm{Ko}}\epsilon^{2/3}k_{n}^{-2/3}$ and
$\nu_{n}=\nu_{*}\sqrt{K_{\mathrm{Ko}}}\epsilon^{1/3}k_{n}^{-4/3}$, where
$k_{n},u_{n}$ are the wavenumber and velocity of shell $n$;
$\nu_{*},K_{\mathrm{Ko}}$ are RG and Kolmogorov’s constants; and $\epsilon$ is
the energy dissipation rate. We find that $\nu_{*}\approx 0.5$ and
$K_{\mathrm{Ko}}\approx 1.7$, consistent with earlier RG works on Navier-
Stokes equation. We verify the theoretical predictions using numerical
simulations.

##  I Introduction

Renormalization group (RG) analysis has been employed to model turbulence.
Orszag (1973) and Forster  _et al._ (1977) performed one of the first
perturbative renormalization analysis. Yakhot and Orszag (1986) performed
detailed analysis using $\epsilon$ expansion. The other perturbative RG works
are by Zhou  _et al._ (1989); Zhou (2010), McComb and Shanmugasundaram (1983,
1985); McComb (2004, 2014), Eyink (1993), Martin  _et al._ (1973),
Bhattacharjee (1991), and Adzhemyan  _et al._ (1999). Among these works,
McComb, Zhou, and coworkers employed self-consistent RG (using “dressed
Green’s function”) that has nonperturbative features. The above set of works
show that the renormalized turbulent viscosity $\nu(k)\sim k^{-4/3}$, where
$k$ is the wavenumber.

Recently, researchers have employed exact renormalization group equation
(ERGE) to turbulence Wilson (1971, 1983); Polchinski (1984). Here, either
sharp or smooth filter is employed during coarsening. A more formal
implementation of ERGE is via functional renormalization group (FRG).
Tomassini (1997), Fontaine  _et al._ (2022), and Canet (2022) employed FRG to
hydrodynamic turbulence and shell model. They derived formulas for the
velocity correlations and multiscaling exponents. For Navier-Stokes equation,
Canet (2022) reported $\nu(k)\approx k^{-1}$, rather than $\nu(k)\sim
k^{-4/3}$. Fedorenko  _et al._ (2013) performed FRG to decaying Burgers,
hydrodynamic, and quasi-geostrophic turbulence. Among many results, Fedorenko
_et al._ (2013) showed that for hydrodynamic turbulence, the second-order
structure function scales as $l$ (the distance between two points), rather
than Kolmogorov’s predictions of $l^{2/3}$.

Mejía-Monasterio and Muratore-Ginanneschi (2012) performed nonperturbative
renormalization group analysis of stochastic Navier-Stokes equation with
power-law forcing. Here, they renormalized the viscosity, the forcing
amplitude, and the coupling constants. Using field-theoretic tools, Biferale
_et al._ (2017) constructed optimal subgrid closure for the shell models; they
related the closure scheme to  large-eddy simulations. In addition, Eyink
(1993) used _operator product expansion_ (OPE) and discovered multiscaling for
the shell model. Some other notable field-theoretic works (not RG) on
turbulence are Kraichnan (1959); Leslie (1973); Fairhall  _et al._ (1997);
Falkovich  _et al._ (2001).

In this paper, we employ RG scheme based on the differential equation, as in
Yakhot and Orszag (1986); McComb (2004); Zhou  _et al._ (1989). Note that the
shell model involves discrete wavenumbers, hence its renormalization does not
involve complex integration, as in hydrodynamic turbulence. For inviscid shell
model, our RG procedure yields $\nu(k)=0$ as the solution of the RG equation,
which is similar to the Gaussian fixed point of Wilson $\phi^{4}$ theory
Wilson and Kogut (1974). We verify several RG predictions using numerical
simulation of the shell model. We use temporal autocorrelation function for
the velocity field to compute the renormalized viscosity Sanada and
Shanmugasundaram (1992); Verma  _et al._ (2020a).

In one of the important works on hydrodynamic turbulence, Kraichnan (1964)
argued that large-scale structures sweep the small-scale fluctuations; this
phenomenon is referred to as sweeping effect. These interactions are naturally
multiscale (across many wavenumbers). Note, however, that multiscale
interactions are absent in the shell models, which has local interactions
among the wavenumber shells. Hence, we expect that sweeping effect may be
suppressed in the shell model. This is precisely what we observe in our RG
calculation of the shell model.

In this paper, we compute the renormalized viscosity in the shell model using
momentum-space RG proposed by Wilson Wilson and Kogut (1974) (see Sec. II).
Here, we assume that the coarse-grained velocity field is random satisfying
time-stationarity. Our calculation does not require quasi-Gaussian
approximation for the velocity field. In Sec. III, we compute the energy flux
of the shell model; here, we assume the velocity field to be quasi-Gaussian.
The flux calculation enables us to compute Kolmogorov’s constant.
Interestingly, our predictions for the shell model are quite close to those
for the Navier-Stokes equation. In Sec. IV, we extend our RG calculation to
show that sweeping effect is suppressed in the shell model.

In Sec. V, we describe how we verify the theoretical predictions using
numerical simulations. We observe that the numerical results are in good
agreement with the theoretical predictions. In Sec. VI, we compare our results
with those from earlier works. We conclude in Sec. VII.

##  II Renormalization of viscosity

The Sabra shell model is L’vov  _et al._ (1998, 1999); Biferale (2003);
Constantin  _et al._ (2006); Ditlevsen (2010); Plunian  _et al._ (2012)

$\displaystyle\frac{du_{n}}{dt}+\bar{\nu}k_{n}^{2}u_{n}$ $\displaystyle=$
$\displaystyle-i\lambda[a_{1}k_{n}u_{n+1}^{*}u_{n+2}+a_{2}k_{n-1}u_{n-1}^{*}u_{n+1}$
(1) $\displaystyle~{}~{}~{}-a_{3}k_{n-2}u_{n-1}u_{n-2}]+f_{n},$

where $u_{n}$ represents the velocity field for the shell $n$; $\bar{\nu}$ is
the microscopic kinematic viscosity; $a_{1},a_{2},a_{3}$ are constants with
$a_{1}+a_{2}+a_{3}=0$; and $k_{n}=k_{0}b^{n}$ with $b$ as a constant. In this
paper, we choose $a_{1}=1,a_{2}=-1+1/b$, $a_{3}=-1/b$, and $b$ in the range
$(1.2,2)$. Here, $f_{n}$ represents the forcing, which is employed at small
$n$’s. This forcing injects energy at large scales that cascades to small
scales as the energy flux. Note that triadic interactions of hydrodynamic
turbulence are modelled better with Sabra model than GOY model L’vov  _et al._
(1998).

The coupling constant (coefficient of the nonlinear term) $\lambda$ is not
renormalized due to the Galilean invariance Forster  _et al._ (1977); McComb
(1990, 2014), and we set $\lambda=1$. Refer to Appendix A for details. In
addition, we consider $u_{n}$ to be random, as in fully-developed turbulence,
rather than introducing a separate noise term in the inertial range McComb and
Shanmugasundaram (1983); Zhou  _et al._ (1988); McComb (1990); Zhou  _et al._
(1988); Zhou (2010). Thus, we avoid noise renormalization. In this self-
consistent approach, we renormalize only the viscosity.

Figure 1:  Division of wavenumber shells into $<$ and $>$ partitions during
the computation of $\nu_{n}$ at $k=k_{n}$. Under the coarse-graining, the
$u^{<}$ variables are unaltered, whereas $u^{>}$ variables are averaged out.

Following Wilson Wilson (1983), we coarse-grain the system over a wavenumber
shell, and compute the consequent correction to the viscosity. The wavenumber
space is already divided in the shell model of turbulence, which makes the
computation simpler than that for Navier-Stokes equation. The locality of
interactions too simplifies the RG calculation. We denote the renormalized
viscosity at wavenumber $k_{n}$ by $\nu_{n}$.

Renormalization is often performed in $({\bf k},\omega)$ space. However, for
the shell model, the renormalization calculation in $({\bf k},t)$ space is
concise and convenient. Hence, we adopt this scheme. In Appendix B, we will
briefly discuss the renormalization of the shell model in $({\bf k},\omega)$
space.

For computing the renormalized viscosity at $k_{n}$ in the inertial range
where $f_{n}=0$, we coarse-grain the system by averaging over
$u_{n+1}^{*>}(t)$ and $u^{>}_{n+2}(t)$ (see Fig. 1). Following RG convention,
we label the variables to be averaged using $>$ symbol, whereas those to be
retained using $<$ symbol. Under this notation,

$\displaystyle\left(\frac{d}{dt}+\bar{\nu}k_{n}^{2}\right)u^{<}_{n}(t)$
$\displaystyle=$ $\displaystyle-i[a_{1}k_{n}u_{n+1}^{*>}(t)u^{>}_{n+2}(t)$ (2)
$\displaystyle+a_{2}k_{n-1}u_{n-1}^{*>}(t)u^{>}_{n+1}(t)$ $\displaystyle-
a_{3}k_{n-2}u^{<}_{n-1}(t)u^{<}_{n-2}(t)].$

The variables with $<$ superscript remain unaltered under coarse-graining.
However, $u^{>}_{n+1}(t)$ and $u^{>}_{n+2}(t)$ variables are assumed to be
random with zero mean. Note that $u^{>}$ variables need not be Gaussian. Under
these assumptions,

$\displaystyle\left\langle u^{*<}_{n-1}(t)u^{>}_{n+1}(t)\right\rangle$
$\displaystyle=$ $\displaystyle u^{*<}_{n-1}(t)\left\langle
u^{>}_{n+1}(t)\right\rangle=0,$ (3) $\displaystyle\left\langle
u^{<}_{n-2}(t)u^{<}_{n-1}(t)\right\rangle$ $\displaystyle=$ $\displaystyle
u^{<}_{n-2}(t)u^{<}_{n-1}(t).$ (4)

Based on the above simplification,

$\displaystyle\left(\frac{d}{dt}+\bar{\nu}k_{n}^{2}\right)u^{<}_{n}(t)-ia_{3}k_{n-2}u^{<}_{n-1}(t)u^{<}_{n-2}(t)$
$\displaystyle=-ia_{1}k_{n}\left\langle
u_{n+1}^{*>}(t)u^{>}_{n+2}(t)\right\rangle.$ (5)

To compute the right-hand side (RHS) of Eq. (5), we evaluate $u_{n+1}^{*>}(t)$
and $u^{>}_{n+2}(t)$ using Green’s function technique. For example,

$\displaystyle u_{n+2}^{>}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{t}dt^{\prime}G_{n+2}(t-t^{\prime})\times(-i)[a_{1}k_{n}u_{n+3}^{*>}(t^{\prime})u^{>}_{n+4}(t^{\prime})$
(6)
$\displaystyle+a_{2}k_{n-1}u_{n+1}^{*>}(t^{\prime})u^{>}_{n+3}(t^{\prime})$
$\displaystyle-a_{3}k_{n}u^{<}_{n}(t^{\prime})u^{>}_{n+1}(t^{\prime})],$

where $G_{n+2}(t-t^{\prime})$ is the Green’s function. Note, however, that
$u_{n+3}^{*>}(t)$ and $u^{>}_{n+4}(t)$ are absent at this stage. Hence,

$\displaystyle u_{n+2}^{>}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{t}dt^{\prime}G_{n+2}(t-t^{\prime})ia_{3}k_{n}u^{<}_{n}(t^{\prime})u^{>}_{n+1}(t^{\prime}).$
(7)

Substitution of Eq. (7) in the RHS of Eq. (5) yields

$\displaystyle I_{1}$ $\displaystyle=$
$\displaystyle\int_{0}^{t}dt^{\prime}G_{n+2}(t-t^{\prime})a_{1}a_{3}k_{n}^{2}u^{<}_{n}(t^{\prime})\left\langle
u^{*>}_{n+1}(t)u^{>}_{n+1}(t^{\prime})\right\rangle$ (8) $\displaystyle=$
$\displaystyle
a_{1}a_{3}k_{n}^{2}\int_{0}^{t}dt^{\prime}G_{n+2}(t-t^{\prime})\bar{C}_{n+1}(t-t^{\prime})u^{<}_{n}(t^{\prime}),$

where $\bar{C}_{n+1}(t-t^{\prime})$ is unequal time correlation.

In self-consistent RG procedure, it is assumed that the decay rates of Green’s
and correlation functions are determined by the renormalized viscosity Leslie
(1973); McComb (2014). Hence,

$\displaystyle G_{n}(t-t^{\prime})$ $\displaystyle=$
$\displaystyle\theta(t-t^{\prime})\exp[-\nu_{n}k_{n}^{2}(t-t^{\prime})],$ (9)
$\displaystyle\bar{C}_{n}(t-t^{\prime})$ $\displaystyle=$ $\displaystyle
C_{n}(t)\exp[-\nu_{n}k_{n}^{2}(t-t^{\prime})],$ (10)

where $C_{n}(t)$ is equal-time correlation ($t=t^{\prime}$), and
$\theta(t-t^{\prime})$ is the step function. Note that $G_{n}(\tau)$ and
$\bar{C}_{n}(\tau)$ decay with a time scale of
$\tau_{c}=(\nu_{n}k_{n}^{2})^{-1}$. Equations (9, 10) are valid for
$\tau<\tau_{c}$, after which $C_{n}$ and $G_{n}$ decay rapidly to zero Orszag
(1973); Pope (2000); Sanada and Shanmugasundaram (1992); Verma  _et al._
(2020a); Zhou (2021); Verma  _et al._ (2020a).

Substitution of $G_{n}(t-t^{\prime})$ and $\bar{C}_{n}(t-t^{\prime})$ of Eqs.
(9, 10) in Eq. (8) yields

$\displaystyle I_{1}$ $\displaystyle=$ $\displaystyle
a_{1}a_{3}k_{n}^{2}C_{n+1}(t)$
$\displaystyle\times\int_{0}^{t}dt^{\prime}\exp[-(\nu_{n+1}k_{n+1}^{2}+\nu_{n+2}k_{n+1}^{2})(t-t^{\prime})]u^{<}_{n}(t^{\prime}),$

Now, we employ Markovian approximation, according to which the integral of Eq.
(II) gets maximal contributions from $t^{\prime}$ near $t$ Orszag (1973). This
is possible when $\nu_{n}k_{n}^{2}\gg 1$ Orszag (1973). Since the integral is
peaked near $t=t^{\prime}$, $u_{n}(t^{\prime})\rightarrow u_{n}(t)$ and

$\displaystyle I_{1}$ $\displaystyle=$
$\displaystyle\frac{a_{1}a_{3}k_{n}^{2}C_{n+1}(t)}{\nu_{n+1}k_{n+1}^{2}+\nu_{n+2}k_{n+2}^{2}}u^{<}_{n}(t).$
(12)

Such assumptions are made in Eddy-damped Quasi-normal Markovian (EDQNM)
approximation of hydrodynamic turbulence Orszag (1973).

The RHS of Eq.(5) has another contribution to $\nu_{n}$, which is computed by
expanding $u^{*>}_{n+1}(t)$ using Green’s function. Following similar approach
as above, we compute new term as

$\displaystyle I_{2}$ $\displaystyle=$
$\displaystyle\frac{a_{1}a_{2}k_{n}^{2}C_{n+2}(t)}{\nu_{n+1}k_{n+1}^{2}+\nu_{n+2}k_{n+2}^{2}}u^{<}_{n}(t).$
(13)

Figure 2: Feynman diagrams associated with the viscosity renormalization.
These diagrams are related to the RHS of Eq. (5).

The Feynman diagrams associated with $I_{1}$ and $I_{2}$ are exhibited in Fig.
2. Here, the loop-diagrams represent the self-energy in which the wavy and
solid lines are the Green’s function and correlation function respectively.

These calculations reveal that the RHS of Eq.(5) is proportional to
$u^{<}_{n}$. Hence, the prefactors of $I_{1}$ and $I_{2}$ will provide
corrections to $\bar{\nu}$ to yield $\nu_{n}$. That is,

$\displaystyle\nu_{n}k_{n}^{2}=\bar{\nu}k_{n}^{2}-\frac{a_{1}k_{n}^{2}[a_{3}C_{n+1}(t)+a_{2}C_{n+2}(t)]}{\nu_{n+1}k_{n+1}^{2}+\nu_{n+2}k_{n+2}^{2}}.$
(14)

Note, however, that $\bar{\nu}\ll\nu_{n}$. Hence,

$\displaystyle\nu_{n}k_{n}^{2}=-\frac{a_{1}k_{n}^{2}[a_{3}C_{n+1}(t)+a_{2}C_{n+2}(t)]}{\nu_{n+1}k_{n+1}^{2}+\nu_{n+2}k_{n+2}^{2}}.$
(15)

Note that we compute renormalized viscosity at the corresponding coarse-
graining step. At the present level, $\nu_{n+1},\nu_{n+2},...$ have been
computed already, whereas, $\nu_{n-1},\nu_{n-2},...$ would be computed at
subsequent stages. Also note that during the computation of $\nu_{n-1}$,
$u^{>}_{n}$ and $u^{>}_{n+1}$ would belong to $>$ shells.

In Eq. (15), $\nu_{n}$ and $C_{n}$ are both unknowns. RG equation for Navier-
Stokes equation too has a similar implicit form. Zhou  _et al._ (1988), and
McComb and Shanmugasundaram (1983) employed self-consistent procedure to solve
such an implicit equation (also see McComb (1990); Zhou (1993, 2010)).
Following these authors, we attempt the following functions for $C_{n}(t)$ and
$\nu_{n}$, which are inspired by Kolmogorov’s theory of turbulence:

$\displaystyle C_{n}(t)$ $\displaystyle=$ $\displaystyle
K_{\mathrm{Ko}}\epsilon^{2/3}k_{n}^{-2/3},$ (16)
$\displaystyle\nu_{n}k_{n}^{2}$ $\displaystyle=$
$\displaystyle\nu_{*}K_{\mathrm{Ko}}^{1/2}\epsilon^{1/3}k_{n}^{2/3},$ (17)

where $K_{\mathrm{Ko}}$ is Kolmogorov’s constant, $\epsilon$ is the viscous
dissipation rate, and $\nu_{*}$ is the RG constant associated with $\nu_{n}$.
Substitution of the above in Eq. (15) yields

$\nu_{*}^{2}=-\frac{a_{1}(a_{3}b^{-2/3}+a_{2}b^{-4/3})}{b^{2/3}+b^{4/3}}.$
(18)

In Fig. 3 we plot $\nu_{*}$ for $b$ ranging from 1.2 to 2.0. Here,
$\nu_{*}\approx 0.5$, in particular, $\nu_{*}\approx 0.48$ for $b=1.5$. The
$\nu_{*}$ computed above are remarkably close to that for Navier-Stokes
equation McComb and Shanmugasundaram (1983); Zhou  _et al._ (1988); McComb
(1990); Zhou (1993, 2010); Verma (2020, 2004), which gives credence to the RG
computation described in this paper.

Figure 3:  For the shell model with various $b$’s, the RG constant $\nu_{*}$
computed using RG (solid blue curve) and using numerical simulations (blue
circles). Also, Kolmogorov’s constant $K_{\mathrm{Ko}}$ computed using field
theory (solid red curve) and using numerical simulations (red squares). The
analytical and numerical $\nu_{*}$’s match quite well, but numerical
$K_{\mathrm{Ko}}$ is around 1.6 times smaller than the analytical counterpart.

It is important to note that the above derivation does not require quasi-
Gaussian assumption for $u^{>}$ variables. We only need to assume time-
stationarity for these variables. In addition, we approximate $\left\langle
u^{<}u^{>}\right\rangle=u^{<}\left\langle u^{>}\right\rangle=0$, rather than
expanding it further. These assumptions and local interactions in the shell
model provide simplification in comparison to the RG calculations for the
Navier-Stokes equation McComb and Shanmugasundaram (1983); Zhou  _et al._
(1988); McComb (1990); Zhou (1993, 2010).

Equation (17) yields

$\nu_{n+1}/\nu_{n}=(k_{n+1}/k_{n})^{-4/3}=b^{-4/3}.$ (19)

As is customary in quantum field theory Peskin and Schroeder (1995), we make a
change of variable as $b=\exp(l)$, with which

$\nu_{n+1}=\nu_{n}\exp(-4l/3)\approx\nu_{n}[1-4l/3],$ (20)

when $b\rightarrow 1$ or $l\rightarrow 0$. Hence,

$\frac{d\nu}{dl}\approx-\frac{4}{3}\nu.$ (21)

Therefore, $\nu_{n}$ increases with the decrease of $k_{n}$, akin to running
coupling constant in quantum chromodynamics. Note, however, that $\nu_{n}$ is
not the coupling constant; instead, it is the coefficient of the viscous term,
which is linear (analogous to mass term in quantum field theory). We remark
that the scaling of Eqs. (16, 17, 21) breaks down when $k\rightarrow 1/L$,
where $L$ is the system size.

The dominant frequency at $k=k_{n}$ is

$\omega_{n}\sim\nu_{n}k_{n}^{2}\sim\epsilon^{1/3}k_{n}^{2/3}.$ (22)

For small $k_{n}$, $\omega_{n}\rightarrow 0$. This is one of the assumptions
of RG schemes in $({\bf k},\omega)$ space. Refer to Appendix B for details.

For $\bar{\nu}=0$, $f_{n}=0$, and $\delta$-correlated (white noise) initial
condition, $u_{n}$ remains $\delta$-correlated, as in Euler turbulence
Kraichnan (1964); Verma  _et al._ (2020b); Verma and Chatterjee (2022).
Therefore, $\left\langle u_{n+1}^{*>}(t)u_{n+2}^{>}(t)\right\rangle=0$ [see
Eq. (5)], leading to no correction or renormalization of the viscosity. Thus,
$\nu_{n}=0$ for the inviscid shell model. This solution corresponds to the
Gaussian fixed point in Wilson’s $\phi^{4}$ theory Wilson and Kogut (1974).

In Sec. III, we will compute the energy flux for the shell model using field-
theoretic techniques.

##  III Energy flux computation

In this section, we compute the energy flux for the shell model
perturbatively. The energy flux at $k=k_{n}$ is defined as Ditlevsen (2010);
Biferale (2003); Verma (2019)

$\displaystyle\Pi_{n}$ $\displaystyle=$ $\displaystyle
2a_{3}k_{n-1}\Im[\left\langle
u^{*}_{n-1}(t)u^{*}_{n}(t)u_{n+1}(t)\right\rangle]$ (23)
$\displaystyle-2a_{1}k_{n}\Im[\left\langle
u^{*}_{n}(t)u^{*}_{n+1}(t)u_{n+2}(t)\right\rangle].$

We compute $\left\langle\Pi_{n}\right\rangle$ by averaging Eq. (23) under the
assumption that $u_{n}(t)$’s in the inertial range are quasi-Gaussian with
zero mean, an assumption used in eddy-damped quasi-normal Markovian (EDQNM)
approximation and in direct interaction approximation (DIA) Kraichnan (1959);
Orszag (1973). To zeroth order, $\left\langle\Pi_{n}\right\rangle=0$, which is
the energy flux for Euler turbulence; this flux corresponds to the Gaussian
fixed point, $\nu=0$.

Figure 4: Feynman diagrams associated with the first term of Eq. (23). Figure
5: Feynman diagrams associated with the second term of Eq. (23).

However, $\left\langle\Pi_{n}\right\rangle\neq 0$ to the first order of
perturbation. The Feynman diagrams associated with the first order in
perturbation for the first and second terms of Eq. (23) are exhibited in Figs.
4 and 5 respectively. Let us analyze the expansion of the first Feynman
diagram of Fig. 4. Here, $u_{n+1}(t)$ has been expanded as

$\displaystyle u_{n+1}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{t}dt^{\prime}G_{n+1}(t-t^{\prime})[-ia_{1}k_{n+1}u_{n+2}^{*}(t^{\prime})u_{n+3}(t^{\prime})$
(24) $\displaystyle\hskip 36.135pt-
ia_{2}k_{n}u_{n}^{*}(t^{\prime})u_{n+2}(t^{\prime})$ $\displaystyle\hskip
36.135pt+ia_{3}k_{n-1}u_{n}(t^{\prime})u_{n-1}(t^{\prime})],$

substitution of which in the first term of Eq. (23), or in the first Feynman
diagram of Fig. 4, yields

$\displaystyle I_{3}$ $\displaystyle=$ $\displaystyle
2a_{3}^{2}k_{n-1}^{2}\int_{0}^{t}dt^{\prime}\exp[-\nu_{n+1}k_{n+1}^{2}(t-t^{\prime})]$
(25) $\displaystyle\times\Im[i\left\langle
u_{n}^{*}(t)u^{*}_{n-1}(t)u_{n}(t^{\prime})u_{n-1}(t^{\prime})\right\rangle]$
$\displaystyle=$ $\displaystyle
2a_{3}^{2}k_{n-1}^{2}\int_{0}^{t}dt^{\prime}\exp[-\nu_{n+1}k_{n+1}^{2}(t-t^{\prime})]$
$\displaystyle\times\left\langle
u_{n}^{*}(t)u_{n}(t^{\prime})\right\rangle\left\langle
u^{*}_{n-1}(t)u_{n-1}(t^{\prime})\right\rangle$ $\displaystyle=$
$\displaystyle
2a_{3}^{2}k_{n-1}^{2}\frac{C_{n}C_{n-1}}{\nu_{n-1}k_{n-1}^{2}+\nu_{n}k_{n}^{2}+\nu_{n+1}k_{n+1}^{2}}.$

In the above derivation, we use the following properties:

  1. 1.

$\left\langle
u_{n}(t)u^{*}_{m}(t^{\prime})\right\rangle=\delta_{n,m}C_{n}(t)\exp[-\nu_{n}k_{n}^{2}(t-t^{\prime})]$.

  2. 2.

$\left\langle abcd\right\rangle=\left\langle ab\right\rangle\left\langle
cd\right\rangle+\left\langle ac\right\rangle\left\langle
bd\right\rangle+\left\langle ad\right\rangle\left\langle bc\right\rangle$ when
$a,b,c,d$ are Gaussian variables.

Using similar analysis, we derive the other integrals of the energy flux as

$\displaystyle I_{4}$ $\displaystyle=$ $\displaystyle
2a_{2}a_{3}k_{n-1}^{2}C_{n-1}C_{n+1}/\mathrm{denr1.},$ (26) $\displaystyle
I_{5}$ $\displaystyle=$ $\displaystyle
2a_{1}a_{3}k_{n-1}^{2}C_{n}C_{n+1}/\mathrm{denr1.},$ (27) $\displaystyle
I_{6}$ $\displaystyle=$
$\displaystyle-2a_{1}^{2}k_{n}^{2}C_{n+1}C_{n+2}/\mathrm{denr2.},$ (28)
$\displaystyle I_{7}$ $\displaystyle=$
$\displaystyle-2a_{1}a_{2}k_{n}^{2}C_{n}C_{n+2}/\mathrm{denr2.},$ (29)
$\displaystyle I_{8}$ $\displaystyle=$
$\displaystyle-2a_{1}a_{3}k_{n}^{2}C_{n}C_{n+1}/\mathrm{denr2.},$ (30)

with

$\displaystyle\mathrm{denr1.}=\nu_{n-1}k_{n-1}^{2}+\nu_{n}k_{n}^{2}+\nu_{n+1}k_{n+1}^{2},$
(31)
$\displaystyle\mathrm{denr2.}=\nu_{n}k_{n}^{2}+\nu_{n+1}k_{n+1}^{2}+\nu_{n+2}k_{n+2}^{2}.$
(32)

The integrals $I_{3},I_{4},I_{5}$ correspond to the first term of Eq. (23),
whereas $I_{6},I_{7},I_{8}$ correspond to the second term of Eq. (23). By
adding $I_{3}$ to $I_{8}$ and using $k_{n}=k_{0}b^{n}$, we derive

$\left\langle\Pi_{n}\right\rangle=\epsilon=\frac{K_{\mathrm{Ko}}^{3/2}}{\nu_{*}}\frac{\mathrm{numr}}{1+b^{2/3}+b^{4/3}},$
(33)

where

$\displaystyle\mathrm{numr}$ $\displaystyle=$ $\displaystyle
2a_{3}b^{-4/3}(a_{1}b^{-2/3}+a_{2}+a_{3}b^{2/3})$ (34)
$\displaystyle-2a_{1}(a_{1}b^{-2}+a_{2}b^{-4/3}+a_{3}b^{-2/3}).$

Equation (33) reveals that the energy flux is independent of wavenumber,
consistent with Kolmogorov’s theory of turbulence Kolmogorov (1941a, b);
Frisch (1995). Using Eq. (33), we compute $K_{\mathrm{Ko}}$ and plot it in
Fig. 3. We observe $K_{\mathrm{Ko}}$ to be a weak function of $b$. In
particular, for $b=1.5$, $K_{\mathrm{Ko}}\approx 1.71$, which is close to the
theoretical, experimental, and numerical values of Kolmogorov’s constant
Frisch (1995); McComb (1990).

##  IV Sweeping Effect in Shell Model

Kraichnan (1964) showed that large-scale flow structures sweep smaller ones, a
phenomenon called sweeping effect. Here, large-scale velocity structures
interact with small-scale ones. Kraichnan (1964) observed that the sweeping
effect leads to $k^{-3/2}$ energy spectrum, rather than usual $k^{-5/3}$
spectrum. To overcome this discrepancy, Kraichnan (1965) proposed
Lagrangian‐History Closure Approximation for Turbulence. Note that the shell
model involves local interactions, thus drastically reduce the sweeping
effect.

To test the sweeping effect in the field-theoretic calculation of shell model,
we introduce a term $iU_{0}k_{n}u_{n}$ in the left-hand side of Eq. (1), where
$U_{0}$, a constant, represents the mean flow. Under renormalization, the
above term appears as $iU_{n}k_{n}u_{n}^{<}$, where $U_{n}$ represents the
renormalized parameter corresponding to $U_{0}$. With $U_{n}$, the RG flow
equation [Eq. (15)] gets transformed to

$\displaystyle
iU_{n}k_{n}+\nu_{n}k_{n}^{2}=-\frac{a_{1}k_{n}^{2}[a_{3}C_{n+1}+a_{2}C_{n+2}]}{\mathrm{denr.}},$
(35)
$\displaystyle\mathrm{denr.}=i[U_{n+1}k_{n+1}+U_{n+2}k_{n+2}]+\nu_{n+1}k_{n+1}^{2}+\nu_{n+2}k_{n+2}^{2}.$
(36)

Using dimensional analysis, we argue that

$U_{n}=U_{*}\epsilon^{1/3}k_{n}^{-1/3},$ (37)

substitution of Eqs. (16,17, 37) in Eq. (35) yields

$(iU_{*}+\nu_{*})^{2}=-\frac{a_{1}[a_{3}b^{-2/3}+a_{2}b^{-4/3}]}{b^{2/3}+b^{4/3}}.$
(38)

The only possible solution of Eq. (38) is

$U_{*}=0,~{}\mathrm{or}~{}U_{n}=0,$ (39)

and $\nu_{*}$ is given by the same formula as Eq. (18). Hence, sweeping effect
is absent in the RG calculation of the shell model, and $\nu(k)$ is
independent of $U_{0}$. However, in Sec. V, we show that the numerical results
deviate from the above prediction.

##  V Numerical Verification

To test the predictions of the above field-theoretic calculations, we solve
Sabra shell model, Eq. (1), numerically. We employ 40 shells, $\nu=10^{-6}$,
$U_{0}=0$, and fourth-order Runge-Kutta (RK4) time marching scheme with
$dt=10^{-5}$. The shell model is forced randomly at shells $n=0$ and 1 so as
to provide a constant energy supply rate; we choose $\epsilon=2$ for all our
runs. To test the dependence of $\nu$ and $K_{\mathrm{Ko}}$ on $b$, we vary
$b$ from 1.2 to 2 in the interval of 0.1. We also perform another simulation
with $U_{0}=0.5$ and $b=1.5$ to test the field-theoretic predictions on the
sweeping effect. We carry out the simulations till 2000 eddy turnover times,
and report the energy spectra and fluxes after the system has reached a steady
state.

As expected, for $U_{0}=0$ and finite $\nu$, in the inertial range, the energy
spectrum $C_{n}\sim k^{-2/3}$ and the energy flux $\Pi_{n}\approx\epsilon=2$
Ditlevsen (2010); Biferale (2003); Plunian  _et al._ (2012). See the red
curves of Fig. 6 for an illustration for $b=1.5$. Using Eq. (16), we compute
the Kolmogorov’s constants for various $b$’s and plot them in Fig. 3. We
observe that for $b=1.5$, $K_{\mathrm{Ko}}=1.05$, which is approximately 1.6
times smaller than the theoretically predicted value of 1.71 (for the shell
model). See Fig. 3 for an illustration. This discrepancy between the numerical
and analytical $K_{\mathrm{Ko}}$ is possibly due to various approximations
employed in the theoretical calculations, an issue that needs a closer
investigation.

Figure 6:  For the numerical simulation of shell model with $b=1.5$: (a) plots
of normalized energy spectra $C_{n}\epsilon^{-2/3}k_{n}^{2/3}$ vs. $k_{n}$ for
$U_{0}=0$ (red curve) and $U_{0}=0.5$ (blue curve). (b) The corresponding
energy fluxes $\Pi_{n}$ are shown using the same color convention. We observe
$C_{n}\sim k_{n}^{-2/3}$ and constant $\Pi_{n}$ in the inertial range. The
insets in (a,b) exhibit $C_{n}$ and $\Pi_{n}$ for the $\nu=0$ case
(equilibrium behaviour).

For the $U_{0}=0.5$ run, we again observe Kolmogorov’s spectrum (apart from a
hump) and constant energy flux (blue curves in Fig. 6). Here,
$K_{\mathrm{Ko}}\approx 0.92$. For the special case with $\bar{\nu}=0$ and
white noise initial condition, numerical simulation yields
$C_{n}\approx\mathrm{constant}$ and nearly zero energy flux, consistent with
the field-theoretic predictions. We illustrate the above energy spectrum and
flux in the insets of Fig. 6.

To validate the renormalized viscosity of Eq. (17), we compute $\nu_{n}$
numerically using the normalized correlation function $R_{n}(\tau)$, which is
defined as

$R_{n}(\tau)=\frac{\bar{C}(\tau)}{C_{n}},$ (40)

where $\bar{C}(\tau)$ and $C_{n}$ are the unequal-time and equal-time
correlations respectively (see Eq. (10)).

We observe that the numerically-computed $R_{n}(\tau)$ is real. For small
$\tau$ and inertial range $k_{n}$’s,
$R_{n}(\tau)\approx\exp[-k_{n}^{2/3}\tau]$, which is consistent with Eqs. (10,
17). As illustrated in Fig. 7, for $b=1.5$ and $U_{0}=0$,

$R_{n}(\tau)\approx 1.03\exp(-0.57k_{n}^{2/3}\tau)$ (41)

when $\nu_{n}k_{n}^{2}\tau\lessapprox 1$. A comparison of Eq. (41) with Eq.
(10) reveals that $\nu_{*}\approx
0.57/(K_{\mathrm{Ko}}^{1/2}\epsilon^{1/3})\approx 0.44$, which is in good
agreement with the RG prediction of 0.48 (see Sec. II). However, we cautiously
remark that the numerical $\nu_{*}$ has significant errors.

Figure 7: For the shell-model simulation with $b=1.5$ and $U_{0}=0$, plots of
$R_{n}(\tau)$ vs. $\nu_{n}k_{n}^{2}\tau$ [Eq. (10)]. The chained straight line
represents the best-fit curve $R_{n}(\tau)=1.03\exp(-0.57k_{n}^{2}\tau)$ in
the interval $\nu_{n}k_{n}^{2}\tau=(0.1,0.75)$.

For $U_{0}=0.5$, we compute $R_{n}(\tau)$ and fit it with
$\exp(-\nu_{n}k_{n}^{2}\tau)$. In Fig. 8 we plot $R_{n}(\tau)$ for shells
$n=5$ to 9. We observe that $R_{n}(\tau)$ for $U_{0}=0.5$ is steeper than the
RG predictions. This is contrary to the RG prediction that the mean flow does
not affect the renormalized viscosity. Clearly, the analytical computation
underestimates the dissipation arising due to nonzero $U_{0}$. This issue
needs a closer examination that will be pursued in future.

Figure 8: For the shell-model simulation with $b=1.5$ and $U_{0}=0.5$, plots
of $R_{n}(\tau)$ vs. $\nu_{n}k_{n}^{2}\tau$. The chained straight line
represents the best-fit curve $R_{n}(\tau)=1.06\exp(-0.51k_{n}^{0.85}\tau)$ in
the interval $\nu_{n}k_{n}^{2}\tau=(0.1,0.6)$ for $n=5$. Note that the shells
8 and 9 do not follow the best-fit curve.

In Sec. VI, we compare our results with earlier RG works on the shell model
and hydrodynamic turbulence.

##  VI Comparison with earlier works

There are only a small number of works on renormalization group analysis of
the shell model. Recently, Fontaine  _et al._ (2022) preformed functional
renormalization group (FRG) analysis of the shell model and computed the
multiscaling exponents. They observed that

$\displaystyle C_{n}\sim k_{n}^{-\alpha_{E}},$ (42)
$\displaystyle\nu_{n}k_{n}^{2}\sim k_{n}^{\alpha_{\nu}},$ (43)

with $\alpha_{E}=0.633\pm.004$ and $\alpha_{\nu}=-0.741\pm 0.01$. Substitution
of the above in renormalization group equation [Eq. (15)] yields

$\alpha_{E}+2\alpha_{\nu}=2.$ (44)

Note that Fontaine  _et al._ (2022)’s $\alpha_{E}$ and $\alpha_{\nu}$ satisfy
Eq. (44) to a good approximation. Fontaine  _et al._ (2022) reported that the
proportionality constant for $C_{2}(\tau=0)$, which is $K_{\mathrm{Ko}}$, is
approximately 1.15.

In a different application of field theory, Eyink (1993) employed operator
product expansion to the shell model and computed various correlations and
structure functions. Note, however, that Fontaine  _et al._ (2022) and Eyink
(1993) do not report the RG constant $\nu_{*}$ in their calculation. We remark
that the multiscaling exponents are related to the fluctuations in the energy
flux, i.e., for $\left\langle\Pi_{n}^{2}\right\rangle$ Sreenivasan (1991); Das
and Bhattacharjee (1994); Verma (2020). The self-consistent calculation
presented in this paper may be extendible to the computation of
$\left\langle\Pi_{n}^{2}\right\rangle$.

It is important to compare our predictions on $\nu_{*}$ and $K_{\mathrm{Ko}}$
with the past works on hydrodynamic turbulence. Yakhot and Orszag (1986)
observed that $\nu_{*}=0.39$ and $K_{\mathrm{Ko}}=1.617$. McComb and
Shanmugasundaram (1985) computed that $\nu_{*}\approx 0.40$ and
$K_{\mathrm{Ko}}\approx 1.8$. Zhou  _et al._ (1988) also reported
$\nu_{*}\approx 0.40$. Our field-theoretic computation of the shell model
yields $\nu_{*}\approx 0.50$ and $K_{\mathrm{Ko}}=1.7$, with minor variations
depending on the value of $b$. Using ERG, Tomassini (1997) showed that
$E(k)\sim k^{-1.666\pm 0.001}$, whereas $K_{\mathrm{Ko}}$ lies in the range of
1.124 to 1.785 depending on the chosen function. Clearly, the shell model
predictions are reasonably close to the earlier works on hydrodynamic
turbulence.

There are subtle differences between the RG schemes for the shell model and
hydrodynamic turbulence. The RG procedure for the shell model does not involve
any integral, and hence is simpler than that for HD turbulence. In addition,
we make fewer assumptions in the RG implementation of the shell model. For
example, $u^{>}$ variables are assumed to be time-stationary, but not
necessarily quasi-Gaussian. Note that many past RG works assume that $u^{>}$
is quasi-Gaussian (see e.g., McComb (2014)). In addition, the RG computation
of the shell model is nearly exact. In Eq. (5), we substitute the expansion of
$u^{*>}_{n+1}$ and $u^{>}_{n+2}$ one after the other, and then solve for the
$\nu_{n}$ under Markovian approximation. Also, note that the local
interactions in the shell model suppresses the _sweeping effect_ proposed by
Kraichnan (1964).

We conclude in the next section.

##  VII Conclusions

In this paper, we employ RG analysis to the shell model of turbulence, and
show that a combination of Kolmogorov’s spectrum
$C_{n}=K_{\mathrm{Ko}}\epsilon^{2/3}k_{n}^{-2/3}$ and
$\nu_{n}=\nu_{*}K_{\mathrm{Ko}}^{1/2}\epsilon^{1/3}k_{n}^{-4/3}$ is a solution
of the RG flow equation. Our calculations predict that for $b=1.5$,
$\nu_{*}\approx 0.48$ and $K_{\mathrm{Ko}}\approx 1.71$, which are in good
agreement with the numerical results, except that the numerical
$K_{\mathrm{Ko}}$ is around 1.6 times smaller than the theoretical prediction.
Note that the field-theoretic predictions for the shell model and the Navier-
Stokes equation are close to each other McComb and Shanmugasundaram (1983);
Zhou  _et al._ (1988); McComb (1990); Zhou  _et al._ (1988).

The computation employed in this paper can be easily generalized to the shell
models for scalar and magnetohydrodynamic turbulence. We also believe that the
fluctuations in the energy flux for the shell model could be computed using
the method outlined in this paper.

###### Acknowledgements.

We thank Soumyadeep Chatterjee in help in drawing the Feynman diagrams of the
paper. We thank anonymous referees for useful comments and suggestions. This
work is supported by the project 6104-1 from the Indo-French Centre for the
Promotion of Advanced Research (IFCPAR/CEFIPRA), and the project
PHY/DST/2020455 by Department of Science and Technology, India.

##  Appendix A Galilean invariance leads to non-renormalizibilty of coupling
constant

It can be easily shown that the coupling constant $\lambda$ of Navier-Stokes
equation (NSE) remain unchanged on renormalization due to Galilean invariance
Forster  _et al._ (1977); McComb (1990, 2014). Here, the derivation is
reproduced in brief.

We write the renormalized Navier-Stokes equation as

$\partial_{t}{\bf u}({\bf x},t)+\lambda{\bf u}({\bf x},t)\cdot\nabla{\bf
u}({\bf x},t)=-\nabla p({\bf x},t)+\nu\nabla^{2}{\bf u}({\bf x},t),$ (45)

where ${\bf u}({\bf x},t),p({\bf x},t)$ are the velocity and pressure fields
respectively, $\lambda$ is a measure of the nonlinear interaction, and $\nu$
is the kinematic viscosity. Note that $\lambda=1$ for the original NSE, but it
may get renormalized under scaling.

We consider two reference frames: Lab references frame, where the fluid has
mean velocity ${\bf U}_{0}=U_{0}\hat{x}$, and the moving reference frame,
where the velocity is ${\bf u^{\prime}}({\bf x^{\prime}},t^{\prime})$ with
zero mean. We denote the variables in the lab frames using unprimed variables,
but those in the moving frame using primed variables. The variables in the two
reference frames are related to each other via Galilean transformation, which
is

$x=x^{\prime}+U_{0}t^{\prime};~{}~{}~{}y=y^{\prime};~{}~{}~{}z=z^{\prime};~{}~{}~{}t=t^{\prime};$
(46)
$\partial_{x}=\partial_{x^{\prime}};~{}~{}~{}\partial_{y}=\partial_{y^{\prime}};~{}~{}~{}\partial_{z}=\partial_{z^{\prime}};~{}~{}~{}\partial_{t}=\partial_{t^{\prime}}-U_{0}\partial_{x}^{\prime};$
(47) ${\bf u}({\bf x},t)=U_{0}\hat{x}+{\bf u^{\prime}}({\bf
x^{\prime}},t^{\prime});~{}~{}~{}p({\bf x},t)=p^{\prime}({\bf
x^{\prime}},t^{\prime}),$ (48)

substitution of which in Eq. (45) yields

$\displaystyle\partial_{t^{\prime}}{\bf u^{\prime}}({\bf
x^{\prime}},t^{\prime})+\lambda[U_{0}\hat{x}+{\bf u^{\prime}}({\bf
x^{\prime}},t^{\prime})]\cdot\nabla^{\prime}{\bf u^{\prime}}({\bf
x^{\prime}},t^{\prime})$ $\displaystyle-U_{0}\partial_{x^{\prime}}{\bf
u^{\prime}}({\bf x^{\prime}},t^{\prime})=-\nabla^{\prime}p^{\prime}({\bf
x},t)+\nu\nabla^{\prime 2}{\bf u^{\prime}}({\bf x^{\prime}},t^{\prime}).$ (49)

Note that Eq. (49) is transformed to Eq. (45) in primed variables only if

$\lambda=1.$ (50)

Thus, it has been shown that $\lambda$ is unchanged under RG due to Galilean
invariance. For further discussion, refer to Forster  _et al._ (1977) and
McComb (1990, 2014).

The shell model is written in Fourier space. Hence, it is not possible to
extend the above derivation to the shell model. However, using the analogy
between the shell model and Navier-Stokes equation, it is reasonable to assume
that $\lambda=1$ for the shell model, and that $\lambda$ remains unaltered
under RG operation.

##  Appendix B  Renormalization of shell model in $({\bf k},\omega)$ space

In this Appendix, we will briefly discuss renormalization of the shell model
in $({\bf k},\omega)$ space. Note that the shell model is already divided in
${\bf k}$ space. The forward and inverse Fourier transforms of $u_{n}$ are
defined as follows:

$\displaystyle u_{n}(t)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}u_{n}(\omega)\exp[-i\omega
t],$ (51) $\displaystyle u_{n}(\omega)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}dt\hskip 2.84544ptu_{n}(t)\exp[i\omega
t].$ (52)

Fourier transform of Eq. (1) yields the following equation for
$u^{<}_{n}(\omega)$:

$\displaystyle(-i\omega+\bar{\nu}k_{n}^{2})u_{n}^{<}(\omega)$ $\displaystyle=$
$\displaystyle-i\int\frac{d\omega^{\prime}}{2\pi}[a_{1}k_{n}u_{n+1}^{*>}(\omega^{\prime}-\omega)u^{>}_{n+2}(\omega^{\prime})$
(53)
$\displaystyle+a_{2}k_{n-1}u^{*<}_{n-1}(\omega^{\prime}-\omega)u^{>}_{n+1}(\omega^{\prime})$
$\displaystyle-
a_{3}k_{n-2}u^{<}_{n-2}(\omega-\omega^{\prime})u^{<}_{n-1}(\omega^{\prime})].$

We perform ensemble averaging over $u_{n+1}^{>}$ and $u_{n+2}^{>}$ variables.
Following the method of Sec. II, we arrive at

$\displaystyle\left\langle
u^{*<}_{n-1}(\omega^{\prime}-\omega)u^{>}_{n+1}(\omega^{\prime})\right\rangle$
$\displaystyle=$ $\displaystyle 0,$ (54) $\displaystyle\left\langle
u^{<}_{n-2}(\omega-\omega^{\prime})u^{<}_{n-1}(\omega^{\prime})\right\rangle$
$\displaystyle=$ $\displaystyle
u^{<}_{n-2}(\omega-\omega^{\prime})u^{<}_{n-1}(\omega^{\prime}).$

Consequently, only the first term of Eq. (53) yields a nonzero correction to
the viscosity.

The renormalized viscosity receives contributions from the two Feynman
diagrams of Fig. 2. For the first loop diagram, we expand
$u_{n+2}^{>}(\omega^{\prime})$ as follows

$\displaystyle(-i\omega^{\prime}+\nu_{n+2}k_{n+2}^{2})u_{n+2}^{>}(\omega^{\prime})=$
$\displaystyle-i\int\frac{d\omega^{\prime\prime}}{2\pi}[a_{1}k_{n+2}u_{n+3}^{*>}(\omega^{\prime\prime}-\omega^{\prime})u^{>}_{n+4}(\omega^{\prime\prime})$
$\displaystyle+a_{2}k_{n+1}a_{2}u^{*>}_{n+1}(\omega^{\prime\prime}-\omega^{\prime})u^{>}_{n+3}(\omega^{\prime\prime})$
$\displaystyle-
a_{3}k_{n}u^{<}_{n}(\omega^{\prime}-\omega^{\prime\prime})u^{>}_{n+1}(\omega^{\prime\prime})].$
(56)

Note, however, that $u_{n+3}^{*>}(\omega)$ and $u_{n+4}^{*>}(\omega)$ are
absent at this stage of expansion. Hence, only the last term of Eq. (56)
survives. Therefore,

$\displaystyle
u_{n+2}^{>}(\omega^{\prime})=ia_{3}k_{n}\int\frac{d\omega^{\prime\prime}}{2\pi}\frac{u^{<}_{n}(\omega^{\prime}-\omega^{\prime\prime})u^{>}_{n+1}(\omega^{\prime\prime})}{-i\omega^{\prime}+\nu_{n+2}k_{n+2}^{2}},$
(57)

substitution of which in the RHS of Eq. (53) yields

$\displaystyle X$ $\displaystyle=$ $\displaystyle-
ia_{1}k_{n}\int\frac{d\omega^{\prime}}{2\pi}\left\langle
u_{n+1}^{*>}(\omega^{\prime}-\omega)u^{>}_{n+2}(\omega^{\prime})\right\rangle$
(58) $\displaystyle=$ $\displaystyle
a_{1}a_{3}k_{n}^{2}\int\frac{d\omega^{\prime}}{2\pi}\frac{d\omega^{\prime\prime}}{2\pi}\frac{\left\langle
u^{*>}_{n+1}(\omega^{\prime}-\omega)u^{>}_{n+1}(\omega^{\prime\prime})\right\rangle}{-i\omega^{\prime}+\nu_{n+2}k_{n+2}^{2}}$
$\displaystyle\hskip 72.26999pt\times
u^{<}_{n}(\omega^{\prime}-\omega^{\prime\prime}),$ $\displaystyle=$
$\displaystyle\left[a_{1}a_{3}k_{n}^{2}\int\frac{d\omega^{\prime}}{2\pi}\frac{\left\langle|u^{>}_{n+1}(\omega^{\prime}-\omega)|^{2}\right\rangle}{-i\omega^{\prime}+\nu_{n+2}k_{n+2}^{2}}\right]u^{<}_{n}(\omega).$

Since $\nu_{n}$ is computed for a long time limit, we set $\omega\rightarrow
0$ in the above integral. Hence, the square-bracketed term of Eq. (58) is

$\displaystyle I_{1}$ $\displaystyle=$ $\displaystyle
a_{1}a_{3}k_{n}^{2}\int\frac{d\omega^{\prime}}{2\pi}\frac{\left\langle|u^{>}_{n+1}(\omega^{\prime})|^{2}\right\rangle}{-i\omega^{\prime}+\nu_{n+2}k_{n+2}^{2}}.$
(59)

Now, we employ Wiener-Khinchin theorem to simplify the frequency spectrum as

$\left\langle|u^{>}_{n+1}(\omega^{\prime})|^{2}\right\rangle=\int_{-\infty}^{\infty}d\tau
C_{n+1}(\tau)\exp(-i\omega^{\prime}\tau),$ (60)

where $C_{n+1}(\tau)$ is the correlation function defined in Eq. (10). With
this,

$\displaystyle I_{1}$ $\displaystyle=$ $\displaystyle
a_{1}a_{3}k_{n}^{2}\int_{-\infty}^{\infty}d\tau\bar{C}_{n+1}(\tau)\int\frac{d\omega^{\prime}}{2\pi}\frac{\exp(-i\omega^{\prime}\tau)}{-i\omega^{\prime}+\nu_{n+2}k_{n+2}^{2}}.$

An application of contour integral over the lower part of $\omega^{\prime}$
plane yields

$\displaystyle I_{1}$ $\displaystyle=$ $\displaystyle
a_{1}a_{3}k_{n}^{2}C_{n+1}$ (62)
$\displaystyle\times\int_{0}^{\infty}d\tau\exp[-(\nu_{n+1}k_{n+1}^{2}+\nu_{n+2}k_{n+2}^{2})\tau]$
$\displaystyle=$
$\displaystyle\frac{a_{1}a_{3}k_{n}^{2}C_{n+1}}{\nu_{n+1}k_{n+1}^{2}+\nu_{n+2}k_{n+2}^{2}}.$

The second Feynman diagram of Fig. 2 yields

$\displaystyle
I_{2}=\frac{a_{1}a_{2}k_{n}^{2}C_{n+2}}{\nu_{n+1}k_{n+1}^{2}+\nu_{n+2}k_{n+2}^{2}}.$
(63)

The steps beyond this point are same as those described in Sec. II.

## References

  * Orszag (1973) S. A. Orszag, in  _Les Houches Summer School of Theoretical Physics_ , edited by R. Balian and J. L. Peube (1973) p. 235.
  * Forster  _et al._ (1977) D. Forster, D. R. Nelson, and M. J. Stephen, Phys. Rev. A 16, 732 (1977).
  * Yakhot and Orszag (1986) V. Yakhot and S. A. Orszag, J. Sci. Comput. 1, 3 (1986).
  * Zhou  _et al._ (1989) Y. Zhou, G. Vahala, and M. Hossain, Phys. Rev. A 40, 5865 (1989).
  * Zhou (2010) Y. Zhou, Phys. Rep. 488, 1 (2010).
  * McComb and Shanmugasundaram (1983) W. D. McComb and V. Shanmugasundaram, Phys. Rev. A 28, 2588 (1983).
  * McComb and Shanmugasundaram (1985) W. D. McComb and V. Shanmugasundaram, J. Phys. A: Math. Theor. 18, 2191 (1985).
  * McComb (2004) W. D. McComb,  _Renormalization Methods: A Guide For Beginners_ (Oxford University Press, Oxford, 2004).
  * McComb (2014) W. D. McComb,  _Homogeneous, Isotropic Turbulence: Phenomenology, Renormalization and Statistical Closures_ (Oxford University Press, 2014).
  * Eyink (1993) G. L. Eyink, Phys. Rev. E 48, 1823 (1993).
  * Martin  _et al._ (1973) P. C. Martin, E. D. Siggia, and H. A. Rose, Phys. Rev. A 8, 423 (1973).
  * Bhattacharjee (1991) J. K. Bhattacharjee, Phys. Fluids A 3, 879 (1991).
  * Adzhemyan  _et al._ (1999) L. T. Adzhemyan, N. V. Antonov, and A. N. Vasiliev,  _Field Theoretic Renormalization Group in Fully Developed Turbulence_ (CRC Press, Boca Raton, FL, 1999).
  * Wilson (1971) K. G. Wilson, Phys. Rev. D 3, 1818 (1971).
  * Wilson (1983) K. G. Wilson, Rev. Mod. Phys. 55, 583 (1983).
  * Polchinski (1984) J. Polchinski, Nucl. Phys. B 231, 269 (1984).
  * Tomassini (1997) P. Tomassini, Phys. Lett. B 411, 117 (1997).
  * Fontaine  _et al._ (2022) C. Fontaine, M. Tarpin, F. Bouchet, and L. Canet, arXiv 10.48550/arxiv.2208.00225 (2022).
  * Canet (2022) L. Canet, J. Fluid Mech. 950, P1 (2022).
  * Fedorenko  _et al._ (2013) A. A. Fedorenko, P. L. Doussal, and K. J. Wiese, J. Stat. Mech. Theory Exp. 2013, P04014 (2013).
  * Mejía-Monasterio and Muratore-Ginanneschi (2012) C. Mejía-Monasterio and P. Muratore-Ginanneschi, Phys. Rev. E 86, 016315 (2012).
  * Biferale  _et al._ (2017) L. Biferale, A. A. Mailybaev, and G. Parisi, Phys. Rev. E 95, 043108 (2017).
  * Kraichnan (1959) R. H. Kraichnan, J. Fluid Mech. 5, 497 (1959).
  * Leslie (1973) D. C. Leslie,  _Developments in the theory of turbulence_ (Clarendon Press, Oxford, 1973).
  * Fairhall  _et al._ (1997) A. L. Fairhall, B. Dhruva, V. S. L’vov, I. Procaccia, and K. R. Sreenivasan, Phys. Rev. Lett. 79, 3174 (1997).
  * Falkovich  _et al._ (2001) G. Falkovich, K. Gawȩdzki, and M. Vergassola, Rev. Mod. Phys. 73, 913 (2001).
  * Wilson and Kogut (1974) K. G. Wilson and J. Kogut, Phys. Rep. 12, 75 (1974).
  * Sanada and Shanmugasundaram (1992) T. Sanada and V. Shanmugasundaram, Phys. Fluids A 4, 1245 (1992).
  * Verma  _et al._ (2020a) M. K. Verma, A. Kumar, and A. Gupta, Trans Indian Natl. Acad. Eng. 5, 649 (2020a).
  * Kraichnan (1964) R. H. Kraichnan, Phys. Fluids 7, 1723 (1964).
  * L’vov  _et al._ (1998) V. S. L’vov, E. Podivilov, A. Pomyalov, I. Procaccia, and D. Vandembroucq, Phys. Rev. E 58, 1811 (1998).
  * L’vov  _et al._ (1999) V. S. L’vov, E. Podivilov, and I. Procaccia, EPL 46, 609 (1999).
  * Biferale (2003) L. Biferale, Annu. Rev. Fluid Mech. 35, 441 (2003).
  * Constantin  _et al._ (2006) P. Constantin, B. Levant, and E. S. Titi, Physica D 219, 120 (2006).
  * Ditlevsen (2010) P. D. Ditlevsen,  _Turbulence and Shell Models_ (Cambridge University Press, Cambridge, 2010).
  * Plunian  _et al._ (2012) F. Plunian, R. Stepanov, and P. Frick, Phys. Rep. 523, 1 (2012).
  * McComb (1990) W. D. McComb,  _The physics of fluid turbulence_ (Clarendon Press, Oxford, 1990).
  * Zhou  _et al._ (1988) Y. Zhou, G. Vahala, and M. Hossain, Phys. Rev. A 37, 2590 (1988).
  * Pope (2000) S. B. Pope,  _Turbulent Flows_ (Cambridge University Press, Cambridge, 2000).
  * Zhou (2021) Y. Zhou, Phys. Rep. 935, 1 (2021).
  * Zhou (1993) Y. Zhou, Phys. Fluids 5, 1092 (1993).
  * Verma (2020) M. K. Verma, arXiv:nlin/0510069v2 (2020).
  * Verma (2004) M. K. Verma, Phys. Rep. 401, 229 (2004).
  * Peskin and Schroeder (1995) M. E. Peskin and D. V. Schroeder,  _An Introduction To Quantum Field Theory_ (The Perseus Books Group, Reading, MA, 1995).
  * Verma  _et al._ (2020b) M. K. Verma, S. Bhattacharya, and S. Bhattacharya, arXiv , arXiv:2004.09053 (2020b).
  * Verma and Chatterjee (2022) M. K. Verma and S. Chatterjee, Phys. Rev. Fluids 7, 114608 (2022).
  * Verma (2019) M. K. Verma,  _Energy transfers in Fluid Flows: Multiscale and Spectral Perspectives_ (Cambridge University Press, Cambridge, 2019).
  * Kolmogorov (1941a) A. N. Kolmogorov, Dokl Acad Nauk SSSR 32, 16 (1941a).
  * Kolmogorov (1941b) A. N. Kolmogorov, Dokl Acad Nauk SSSR 30, 301 (1941b).
  * Frisch (1995) U. Frisch,  _Turbulence: The Legacy of A. N. Kolmogorov_ (Cambridge University Press, Cambridge, 1995).
  * Kraichnan (1965) R. H. Kraichnan, Phys. Fluids 8, 575 (1965).
  * Sreenivasan (1991) K. R. Sreenivasan, Annu. Rev. Fluid Mech. 23, 539 (1991).
  * Das and Bhattacharjee (1994) A. Das and J. K. Bhattacharjee, EPL 26, 527 (1994).