# Wall Ion Loss Reduction by Acceleration Zone Shifting in Anode-Layer Hall
Thruster

Rei Kawashima 111Assistant Professor, Department of Aeronautics and
Astronautics<EMAIL_ADDRESS>Member AIAA. The University of
Tokyo, Tokyo, 113-8656, Japan  Yushi Hamada222Graduate Student, Department of
Advanced Energy. Student Member AIAA. The University of Tokyo, Tokyo,
113-8656, Japan  Shu Kawabata333Graduate Student, Department of Aeronautics
and Astronautics. The University of Tokyo, Tokyo, 113-8656, Japan  Kimiya
Komurasaki444Professor, Department of Aeronautics and Astronautics. Member
AIAA. The University of Tokyo, Tokyo, 113-8656, Japan  Hiroyuki
Koizumi555Associate Professor, Department of Advanced Energy. Member AIAA. The
University of Tokyo, Kashiwa, 277-8561, Japan

## Nomenclature

$B_{r}$ = | radial magnetic flux density, T
---|---
$E_{\rm i,w}$ = | ion energy colliding with wall, eV
$E_{\rm ion}$ = | ionization potential, eV
$E_{\rm work}$ = | work function of wall, eV
$e$ = | elementary charge, C
$I_{\rm d}$ = | discharge current, A
$I_{\rm g}$ = | guard ring current, A
$i_{\rm g}$ = | guard ring current density, A/m2
$i_{\rm i,w}$ = | ion wall current density, A/m2
$M$ = | atomic mass of wall material, kg
$P_{\rm i,w}$ = | ion wall energy loss, W
$S$ = | ion collecting surface area, m2
$V_{\rm d}$ = | discharge voltage, V
$V_{\rm g}$ = | guard ring voltage, V
$Y$ = | sputter yield, atoms/ion
$Y_{\rm E}$ = | energy factor of sputter yield, atoms/ion
$Y_{\rm v}$ = | volume sputter yield, mm3/C
$Y_{\theta}$ = | angle factor of sputter yield, –
$\gamma_{\rm e}$ = | electron current factor, –
$\gamma_{\rm SEE}$ = | secondary electron emission yield, –
$\varepsilon$ = | wall erosion rate, µm/h
$\rho$ = | density of wall material, kg/m3

##  1 Introduction

Hall thrusters are electric propulsion systems used in spacecraft. They are
equipped with an axisymmetric discharge channel with a crossed-field
configuration, where radial magnetic and axial electric fields are applied.
Hall thrusters are available in two primary configuration types: a stationary
plasma thruster (SPT) and a thruster with anode layer (TAL). The
distinguishing features of these configurations are discussed in the
literature [1, 2]. Compared with the SPT, the TAL demonstrates high-density
plasma generation and a more compact thruster body with a short discharge
channel. Additionally, the TAL contains metallic channel walls biased at the
cathode potential to repel the electrons approaching the walls, resulting in
reduced electron energy losses to the channel walls. The TAL has the potential
for efficient plasma extraction; however, its superiority has not yet been
demonstrated [3].

The optimal design and operation regime for TAL-type Hall thrusters with xenon
and alternative propellants have been investigated [4, 5]. Hamada et al.
developed a 5-kW TAL characterized by a magnetic field geometry where the peak
of the magnetic flux density is located downstream of the channel exit [6].
This thruster attained a maximum thrust efficiency of 0.60, which is
comparable to the efficiencies of SPTs with the same input power level. A 2-kW
TAL was designed to incorporate discharge channel and magnetic field
geometries similar to those of the 5-kW thruster [7]. The thruster size was
reduced from that of the 5-kW thruster to maintain the flux density of the
propellant xenon gas. Probe diagnostics showed that the bulk plasma properties
of the 5-kW and 2-kW TALs are similar. Furthermore, these TAL thrusters
exhibit primary acceleration zones outside the channel exit, owing to the
abovementioned magnetic field design. This plasma property is different from
that of conventional TALs, in which the acceleration zone exists in the thin
“anode layer” in front of the anode.

In a TAL with the acceleration zone shifted downstream, the plasma–wall
interaction is critical for the energy efficiency and thruster lifetime. If
the potential of the bulk plasma inside the discharge channel is maintained
near the anode potential, a fraction of the ions is inevitably attracted to
the conducting channel walls. The wall ion loss implies a loss of energy for
ionization and acceleration. In addition, the ions colliding with the wall
display high energy because of the large potential drop between the bulk
plasma and wall, leading to wall erosion due to ion bombardment. Wall erosion
is a significant concern in Hall thruster development because it is one of the
primary limiting factors of thruster lifetime. Several technologies have been
proposed for wall erosion reduction in SPTs, including magnetic shielding [8]
and wall-less thrusters [9]. Wall ion loss and erosion must be assessed for
TALs with shifted acceleration zones.

The objective of this study was to measure the wall ion fluxes during TAL
operation with a shifted acceleration zone and evaluate the wall erosion rate.
The ion wall flux was evaluated by measuring the metallic-channel wall
current. A previous study found that the wall current was dependent on the
discharge voltage [7]. In this study, we further investigated the dependence
of the ion wall flux on the magnetic configuration with a detailed analysis of
the guard ring current. A fast estimation method for the wall erosion rate
using the guard ring current is proposed, and the erosion rate of the TAL with
a shifted acceleration zone is evaluated.

##  2 Experimental Setup

The TAL thruster RAIJIN66 was used in the experiment [7]. The thruster
contains a hollow copper anode, and the propellant xenon gas is supplied
through the hollow anode. The metallic channel wall is called the guard ring
and is constructed from stainless steel. A magnetic screen was installed in
the magnetic circuit using recent SPT designs [10]. In the RAIJIN thrusters,
the shifting of the acceleration zone downstream beyond the channel exit was
achieved by using the magnetic screen effectively. Two magnetic field
geometries were tested (with and without the magnetic screen), as illustrated
in Figs. 1 and 2. When the magnetic screen was installed, the magnetic flux
density upstream of the channel exit was trimmed, and thus, the region of
strong magnetic confinement was concentrated downstream. The position of the
peak magnetic flux density was at $z=0$ mm in the case without the magnetic
screen and shifted to $z=1.5$ mm in the case with the magnetic screen. As
displayed in Fig. 1, magnetic flux lines intersect the guard rings in both
cases. Thus, magnetic shielding was not implemented in this thruster.

Figure 1: Magnetic field geometries in cases with and without magnetic screen.
Figure 2: Magnetic flux density distributions along channel centerline.

The circuit diagram during the thruster operation is shown in Fig. 3. In
nominal TAL operation, the guard rings and thruster body, including the front
and back plates, are biased at the cathode potential. The cathode body line
was connected to the ground in this experiment for the stability of thruster
operation. In other words, the cathode-to-ground voltage was set to zero. This
would decrease the potential of the beam relative to the cathode and increase
acceleration potential for ions compared with the case of floating cathode
line. Experimental studies have shown that the cathode-to-ground voltage is
generally negative to the ground, and has a value of 3% to 6% of the discharge
voltage [11]. We assumed that a variation of the acceleration potential over
this range would not change the fundamental characteristics of wall ion
losses. In this experiment, the guard ring was isolated from the other
thruster body parts, and $I_{\rm g}$ was defined as the current flowing into
the inner and outer guard rings. Additionally, a power supply was inserted to
bias the guard rings to a potential lower than that of the cathode. $I_{\rm
d}$ and $I_{\rm g}$ were obtained by measuring the voltage drops across the
0.5-$\Omega$ shunt resistors. Data acquisition was implemented by differential
probes and an oscilloscope.

Figure 3: Circuit diagram of thruster operation and measurement.

To calculate the guard ring current density $i_{\rm g}$, the ion-collecting
surface area $S$ must be defined. A close-up view of the guard ring and
magnetic lines is shown in Fig. 4. The guard rings made from type 304
stainless steel are in contact with the plasma at the surfaces in the r- and
z-directions. The bulk discharge ions generated in the discharge channel flow
into the guard ring from the surface normal to the r-direction. The surface
normal to the z-direction predominantly collect the ions generated in the
charge-exchange collision (CEX) process. It is assumed that the contribution
of the CEX ions to the guard ring current is small compared to that of the
bulk discharge ions. Therefore, the ion-collecting area is defined as the
surface normal to the r-direction, with an axial length of 2.1 mm (distance
between the hollow anode tip and guard ring end) and outer and inner diameters
of 66 mm and 54.1 mm, respectively. Note that neglecting the z-direction
surfaces results in a slight overestimation of $i_{\rm g}$ (and thus of the
wall erosion) because the measured $I_{\rm g}$ includes both the bulk
discharge and CEX ions.

Figure 4: Close-up view of the guard ring.

A vacuum chamber with a diameter of 2 m and length of 3 m was used in the
experiment. During the experiment, the background pressure was maintained
below 6.6$\times$10-3 Pa. The hollow cathode HCN-252 (Veeco Instruments Inc.)
was used as the electron emitter, with a xenon flow rate of 0.29 mg/s. The
background xenon particles inevitably affect the thruster performance and
discharge plasma characteristics. The background pressure effects were
investigated for 1.35-kW SPT-100, and it was reported that the backpressure of
6.7$\times$10-3 Pa caused a non-negligible performance change such as 4%
increase in thrust, because a fraction of the background particles worked as
propellant [12]. In addition, it was reported that the acceleration zone moved
upstream by a few millimeters owing to the backpressure [13]. In this
experiment, it was expected that these background pressure effects would
increase the measured $i_{\rm g}$. Thus, the analysis given in this study may
yield a slightly overestimated wall erosion rate.

The anode mass flow rate of the xenon propellant was set to 3.35 mg/s, and the
discharge voltage was set to either 200 V or 300 V. Under these conditions,
two experiments were carried out: 1) the current–voltage (I-V) characteristics
of the guard rings were obtained to examine how much electron current was
included in the guard ring current, and 2) $i_{\rm g}$ was measured for
different magnetic configurations with various magnetic field strengths. The
systematic errors of the measurement system were estimated to be $\pm 3.7$%
for $I_{\rm d}$ and $I_{\rm g}$, and $\pm 10.2$% for $i_{\rm g}$ based on the
errors in the shunt resistors, differential probes, oscilloscope, and
machining of guard rings.

##  3 Results and Discussion

###  3.1 I-V Characteristics of Guard Ring

The dependence of $I_{\rm d}$ and $I_{\rm g}$ on the guard ring bias voltage
$V_{\rm g}$ is illustrated in Fig. 5. In this experiment, the discharge
voltage and peak magnetic flux density were fixed at 200 V and 23 mT,
respectively. The typical errors for $I_{\rm d}$ and $I_{\rm g}$ are $\pm$0.13
A and $\pm$4.5 mA, respectively. To examine if the guard ring sheath was in
the ion saturation regime, $V_{\rm g}$ was varied over a range of $-50$
V$<V_{\rm g}<$ 0 V. During the sweep of $V_{\rm g}$, $I_{\rm d}$ was
maintained close to 3.34 A, and hence, the effect of the guard ring voltage on
the bulk discharge plasma was insignificant. $I_{\rm g}$ exhibited a slight
dependence on $V_{\rm g}$. Linear fitting was conducted using the data for
$-50$ V$<V_{\rm g}<-20$ V, and the ion saturation current was estimated by
extrapolating the fitting line to $V_{\rm g}$ at 0 V. The dimensionless factor
$\gamma_{\rm e}$ was used to indicate the electron current included in the
guard ring current, which is defined as $1+\gamma_{\rm e}=\left(I_{\rm
g,fit}/I_{\rm g}\right)_{V_{\rm g}=0}$. $\gamma_{\rm e}$ of 0.044 was obtained
from the measurements. The small $\gamma_{\rm e}$ indicates that the sheaths
on the guard rings are in the ion saturation regime and that the guard ring
current can be used to evaluate the wall ion fluxes. A $\gamma_{\rm e}$ value
of 0.044 was used in the wall ion current analysis described in Sec. 3.3.

Figure 5: Discharge and guard-ring currents with varied guard-ring voltages.

###  3.2 Guard Ring Current Measurement Result

Figure 6 displays the guard ring current density corresponding to various
magnetic flux densities for the cases with and without the magnetic screen.
$i_{\rm g}$ is typically dependent on $B_{r}$, indicating that an optimal
magnetic flux density exists in which the guard ring current is minimized. At
$V_{\rm d}=300$ V and $B_{r}=32$ mT, a minimum $i_{\rm g}$ of 12 A/m2 was
obtained during operation with the magnetic screen, whereas a minimum $i_{\rm
g}$ of 98 A/m2 was obtained when the magnetic screen was not used.

Figure 6: Guard ring current densities for various magnetic flux densities.

In TALs, the ions colliding with the guard rings induce secondary electron
emission (SEE). The secondary electrons emitted from the guard ring surface
increase the measured current. The wall ion current density was obtained from
the guard ring current density as $i_{\rm i,w}=\left(1+\gamma_{\rm
e}\right)/\left(1+\gamma_{\rm SEE}\right)\cdot i_{\rm g}$, where $\gamma_{\rm
SEE}$ is the SEE yield. When the acceleration potential of impact ion is in
the range of 10–1000 V, the process of SEE is realized by Auger emission [14].
The ion recombination on the wall surface is a three-body recombination
process where two electrons are involved. When an electron enters the ground
state of the impact ion, the energy corresponding to $E_{\rm ion}$ is
released. Because the first electron is emitted from the surface, the energy
of $E_{\rm work}$ is consumed from the ionization potential. The excess
energy, $E_{\rm ion}-E_{\rm work}$ is transferred to the secondary electron.
The secondary electron must gain an energy greater than $E_{\rm work}$ to move
away from the surface. Thus, the condition for SEE is $E_{\rm ion}-2E_{\rm
work}\geq 0$. In this study, we used an empirical expression to estimate the
SEE yield: $\gamma_{\rm SEE}\approx 0.016\left(E_{\rm ion}-2E_{\rm
work}\right)=0.055$ [14]. Here, $E_{\rm ion}=12.13$ eV is the ionization
potential of xenon, and $E_{\rm work}=4.34$ eV [15] is the work function of
type 304 stainless steel. Using these relationships, the wall ion current
density can be estimated from the guard-ring current density as $i_{\rm
i,w}=0.99i_{\rm g}$. This conversion coefficient was assumed to be constant
for the cases of different magnetic configurations.

The estimated $i_{\rm i,w}$ is 97 A/m2 for the case without the magnetic
screen and 12 A/m2 for the case with the magnetic screen. $i_{\rm i,w}$ is
reduced because the acceleration and ionization zones move downstream when the
magnetic screen is used, and the ions are exhausted without colliding with the
guard rings. Therefore, the wall ion loss in a TAL can be reduced using a
magnetic screen and the optimal magnetic flux density.

###  3.3 Wall Erosion Analysis

The wall erosion rate was estimated based on the measured wall ion current and
sputter yield. This estimation method allows for a quick assessment of the
wall erosion rate in TALs. The sputter yield $Y$ (atoms/ion) was modeled as
$Y=Y_{E}\times Y_{\theta}$, where $Y_{E}$ and $Y_{\theta}$ represent the
functions of the dependences of the incident ion energy and angle,
respectively. Emissive probe measurements of the TAL thruster demonstrated
that the plasma potential was close to the anode voltage inside the discharge
channel [7]. Thus, the acceleration potential of ions colliding with the
cathode-biased wall was assumed to be as high as $V_{\rm d}$. A $Y_{\rm E}$ of
0.55 was assumed for 300-eV xenon ion bombardment on stainless steel walls
based on experimental data [16]. In this analysis, the wall erosion rate was
estimated for the graphite wall material. Carbon graphite exhibits durability
against ion bombardment, and a $Y_{\rm E}$ of 0.08 was assumed for a 300-eV
xenon impact based on published experimental data [17, 18]. $Y_{\theta}$ is a
normalized quantity relative to the case of normal incidence. Several studies
have shown that $Y_{\theta}$ generally increases as the incident angle of the
impact ion increases, and $Y_{\theta}$ has a peak when the incident angle is
around 45–80 deg [19, 20]. In the thruster operation with the magnetic
screens, the axial acceleration inside the channel was negligible compared
with the radial acceleration due to the potential drop in $V_{\rm d}$ toward
the guard rings, concerning the ions colliding with the walls. Therefore, the
normal incidence case was assumed as $Y_{\theta}=1$. In the case without the
magnetic screens, there would be a possibility that the acceleration zone
shifted upstream and the ion incidences were in oblique. However, because the
primary purpose of the present analysis was to analyze the erosion rate of the
thruster with the magnetic screens, we assumed $Y_{\theta}$ of unity for both
magnetic configurations. The volumetric sputter yield is calculated as $Y_{\rm
V}=YM/\left(e\rho\right)$, where $M$ and $\rho$ are the mean atomic weight and
mass density of the wall material. The wall erosion rate can be estimated as
$\varepsilon=Y_{\rm V}i_{\rm i,w}=Y_{E}Y_{\theta}M/\left(e\rho\right)\cdot
i_{\rm i,w}$.

The wall ion energy loss $P_{\rm i,w}$ and erosion rate $\varepsilon$
estimated for RAIJIN66 are listed in Table 2. As a reference, the SPT data in
the unshielded and magnetically shielded configurations are also shown [8].
Here, the power corresponding to the wall ion energy loss was estimated as
$P_{\rm i,w}=i_{\rm i,w}S\left(E_{\rm ion}+E_{\rm i,w}\right)$. In RAIJIN66,
by installing a magnetic screen, the wall ion loss was reduced by
approximately one-eighth. The wall ion energy loss power for the case with the
magnetic screen was reduced to 2.6 W, which was less than 0.3% of the 1 kW
anode input power. The wall erosion rate was also reduced to 1.68 µm/h for the
stainless steel walls. Furthermore, if graphite is used as the guard ring
material, $\varepsilon$ can be expected to be 0.18 µm/h. This value implies
that the wall recession is expected to be 1.8 mm for a 10000 h operation. The
guard rings of RAIJIN66 are 2.0 mm thick, and therefore endure the erosion. It
is concluded that a thruster lifetime exceeding 10,000 h is possible in TALs
with a downstream-shifted acceleration zone.

Table 2:  Energy loss and wall erosion rate estimation for RAIJIN66 and H6 [8]. Thruster |  Wall |  Configuration |  $i_{\rm i,w}$, A m-2 |  $E_{\rm i,w}$, eV |  $Y_{\rm v},\times 10^{-3}$ mm3 C-1 |  $P_{\rm i,w}$, W |  $\varepsilon$, µm/h
---|---|---|---|---|---|---|---
RAIJIN66 |  SUS |  without mag. screen | 97 | 300 | 39 | 22 | 13.7
RAIJIN66 |  SUS |  with mag. screen | 12 | 300 | 39 | 2.6 | 1.68
RAIJIN66 |  C |  with mag. screen | 12 | 300 | 4.3 | 2.6 | 0.18
H6 |  BN |  unshielded | 210 | 143 | 28 | – | 21.1
H6 |  BN |  magnetic shielding | 40 | 36.1 | 0.14 | – | 0.02

Although the wall erosion of RAIJIN66 with a carbon graphite wall was expected
to be minor, the wall erosion rate was greater than that in the magnetically
shielded SPT by one order of magnitude. To further extend the TAL thruster
lifetime, chamfering or curving of the guard ring, such as with the
magnetically shielded thrusters, may be an effective method. In addition, we
identified the optimal operation condition in terms of the wall ion loss and
wall erosion. It was found that this optimal operation condition deviated from
the point of the highest thruster performance. The anode efficiency of
RAIJIN66 was obtained for various operating conditions [7]. For the case with
discharge voltage and anode mass flow rate of 300 V and 3.35 mg/s, $B_{r}$ of
70 mT yielded higher anode efficiency compared with the case of $B_{r}=32$ mT.
This means that there is a trade-off between the wall erosion and anode
efficiency in the present thruster design. Further exploration of the optimal
operating condition and optimization of thruster design are envisaged for a
future study.

##  4 Conclusion

The wall ion loss and erosion rate were investigated for a TAL-type Hall
thruster characterized by an acceleration zone shifted downstream. The guard
rings were electrically separated from the thruster body, and the guard ring
current was measured independently. The I-V characteristics of the guard ring
indicated that the sheath in front of the guard rings was in the ion
saturation regime and that the guard ring current can be used to indicate the
wall ion flux, enabling a fast evaluation of the wall erosion rate. The wall
ion current, wall ion energy loss, and erosion rate were estimated based on
the measured guard ring current.

The results demonstrate that the wall ion loss can be significantly reduced
using a magnetic screen and the optimal magnetic flux density. The wall ion
energy loss was less than 0.3% of the anode input power. The wall erosion rate
was estimated to be 1.68 µm/h and 0.18 µm/h when using stainless steel and
carbon walls, respectively. If carbon material is used, a thruster lifetime
exceeding 10000 h is likely. It has been demonstrated that improved thruster
lifetime can be achieved in the TAL, with the acceleration zone shifted
downstream.

## Funding Sources

This work was supported by JSPS KAKENHI Grant Number JP20H02346.

## References

  * Zhurin et al. [1999] Zhurin, V. V., Kaufman, H. R., and Robinson, R. S., “Physics of Closed Drift Thrusters,” _Plasma Sources Science and Technology_ , Vol. 8, No. 1, 1999, pp. R1–R20.  10.1088/0963-0252/8/1/021.
  * Choueiri [2001] Choueiri, E. Y., “Fundamental Difference between the Two Hall Thruster Variants,” _Physics of Plasmas_ , Vol. 8, No. 11, 2001, pp. 5025–5033.  10.1063/1.1409344.
  * Mazouffre [2016] Mazouffre, S., “Electric Propulsion for Satellites and Spacecraft: Established Technologies and Novel Approaches,” _Plasma Sources Science and Technology_ , Vol. 25, No. 3, 2016, p. 033002.  10.1088/0963-0252/25/3/033002.
  * Yamasaki et al. [2019] Yamasaki, J., Yokota, S., and Shimamura, K., “Performance enhancement of an argon-based propellant in a Hall thruster,” _Vacuum_ , Vol. 167, 2019, pp. 520–523.  10.1016/j.vacuum.2018.09.042.
  * Dukhopelnikov et al. [2021] Dukhopelnikov, D. V., Riazanov, V. A., Shilov, S. O., Manegin, D. S., and Sokolov, R. A., “Investigation of the laboratory model of a thruster with anode layer operating with air and nitrogen-oxygen mixture,” _AIP Conference Proceedings_ , Vol. 2318, No. 1, 2021, p. 040006.  10.1063/5.0036251.
  * Hamada et al. [2017] Hamada, Y., Bak, J., Kawashima, R., Koizumi, H., Komurasaki, K., Yamamoto, N., Egawa, Y., Funaki, I., Iihara, S., Cho, S., Kubota, K., Watanabe, H., Fuchigami, K., Tashiro, Y., Takahata, Y., Kakuma, T., Furukubo, Y., and Tahara, H., “Hall Thruster Development for Japanese Space Propulsion Programs,” _Transactions of the Japan Society for Aeronautical and Space Sciences_ , Vol. 60, No. 5, 2017, pp. 320–326.  10.2322/tjsass.60.320.
  * Hamada et al. [2021] Hamada, Y., Kawashima, R., Bak, J., Komurasaki, K., and Koizumi, H., “Characterization of Acceleration Zone Shifting in an Anode-layer-type Hall Thruster RAIJIN66,” _Vacuum_ , Vol. 186, 2021, p. 110040\.  10.1016/j.vacuum.2020.110040.
  * Hofer et al. [2014] Hofer, R. R., Goebel, D. M., Mikellides, I. G., and Katz, I., “Magnetic Shielding of a Laboratory Hall Thruster. II. Experiments,” _Journal of Applied Physics_ , Vol. 115, No. 4, 2014, p. 043304.  10.1063/1.4862314.
  * Vaudolon et al. [2015] Vaudolon, J., Mazouffre, S., Hénaux, C., Harribey, D., and Rossi, A., “Optimization of a Wall-less Hall Thruster,” _Applied Physics Letters_ , Vol. 107, No. 17, 2015, p. 174103.  10.1063/1.4932196.
  * Reza et al. [2019] Reza, M., Ferrato, E., Faraji, F., Andreussi, T., and Andrenucci, M., “Magnetic Circuit Optimization for Hall Thrusters Design,” 36th International Electric Propulsion Conference, IEPC-2019-A459, 2019.
  * Goebel et al. [2012] Goebel, D. M., Jameson, K. K., and Hofer, R. R., “Hall Thruster Cathode Flow Impact on Coupling Voltage and Cathode Life,” _Journal of Propulsion and Power_ , Vol. 28, No. 2, 2012, pp. 355–363.  10.2514/1.B34275.
  * Diamant et al. [2014] Diamant, K. D., Liang, R., and Corey, R. L., “The Effect of Background Pressure on SPT-100 Hall Thruster Performance,” 50th AIAA/ASME/SAE/ASEE Joint Propulsion Conference, AIAA 2014-3710, 2014.  doi:10.2514/6.2014-3710.
  * MacDonald-Tenenbaum et al. [2019] MacDonald-Tenenbaum, N., Pratt, Q., Nakles, M., Pilgram, N., Holmes, M., and Hargus, W., “Background Pressure Effects on Ion Velocity Distributions in an SPT-100 Hall Thruster,” _Journal of Propulsion and Power_ , Vol. 35, No. 2, 2019, pp. 403–412.  10.2514/1.B37133.
  * Lieberman and Lichtenberg [2005] Lieberman, M. A., and Lichtenberg, A. J., _Principles of Plasma Discharges and Materials Processing_ , 2nd ed., John Wiley & Sons, Hoboken, NJ, 2005.
  * Wilson [1966] Wilson, R. G., “Vacuum Thermionic Work Functions of Polycrystalline Be, Ti, Cr, Fe, Ni, Cu, Pt, and Type 304 Stainless Steel,” _Journal of Applied Physics_ , Vol. 37, No. 6, 1966, pp. 2261–2267.  10.1063/1.1708797.
  * Andersen and Bay [1980] Andersen, H. H., and Bay, H. L., “A Survey of Sputtering-yield Data for Plasma-wall Interaction Calculations,” _Journal of Nuclear Materials_ , Vol. 93-94, 1980, pp. 625–633.  10.1016/0022-3115(80)90184-1.
  * Kenmotsu et al. [2009] Kenmotsu, T., Wada, M., Hyakutake, T., Muramoto, T., and Nishida, M., “Enhanced Sputtering Yields of Carbon Due to Accumulation of Low-energy Xe Ions,” _Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms_ , Vol. 267, No. 8, 2009, pp. 1717–1720.  10.1016/j.nimb.2009.01.127.
  * Boyd and Falk [2001] Boyd, I., and Falk, M., “A Review of Spacecraft Material Sputtering by Hall Thruster Plumes,” 37th Joint Propulsion Conference and Exhibit, AIAA 2001-3353, 2001.  10.2514/6.2001-3353.
  * Yamamura and Shindo [1984] Yamamura, Y., and Shindo, S., “An empirical formula for angular dependence of sputtering yields,” _Radiation Effects_ , Vol. 80, No. 1-2, 1984, pp. 57–72.  10.1080/00337578408222489.
  * Ranjan et al. [2016] Ranjan, M., Sharma, A., Vaid, A., Bhatt, T., Nandalan, V., James, M. G., Revathi, H., and Mukherjee, S., “BN/BNSiO2 sputtering yield shape profiles under stationary plasma thruster operating conditions,” _AIP Advances_ , Vol. 6, No. 9, 2016, p. 095224.  10.1063/1.4964312.


By keeping only the $\order{1/q^{2}}$ contributions to
$\mathcal{D}^{2,1}(x\geq 2,T)$, have found a set of diagrams labelled by:
$\tau$, the length of each of the physical OTOCs; $t_{r}$ for each site
$r\in\\{2,\cdots,x-1\\}$, the positions the ‘$0$’ contraction on site $r$.
Setting $t_{1}=\tau$ and $t_{x}=T-\tau$, we label each diagram by a sequence
$(t_{1},t_{2},\cdots,t_{x-1},t_{x})$. An example diagram, labelled
$(t_{1},t_{2},t_{3},t_{4},t_{5},t_{6})=(5,3,6,2,8,7)$, that contributes to
$\mathcal{D}^{2,1}(x=6,T=12)$ is given below

Figure 4: The contribution
$(t_{1},t_{2},t_{3},t_{4},t_{5},t_{6})=(5,3,6,2,8,7)$ to $\int
dV\mathcal{D}^{(2,1)}(x=6,T=12)$

In fact, for $x>2$, only contributions where $t_{1}\leq
t_{2}<t_{3}<\cdots<t_{x}$ are non-zero. One of the following motifs must
appear in any diagram not satisfying this property, each of which is zero,

$\displaystyle\raisebox{-0.45pt}{\includegraphics[height=56.9055pt]{Diagrams/minus-
square-minus-contraction.png}},$
$\displaystyle\quad\raisebox{-0.45pt}{\includegraphics[height=56.9055pt]{Diagrams/plus-
triangle-plus-contraction.png}},$
$\displaystyle-\raisebox{-0.45pt}{\includegraphics[height=56.9055pt]{Diagrams/square-
square-minus-contraction.png}}.$ (154)

In the non-zero diagrams (i.e., those for which $t_{1}\leq
t_{2}<t_{3}<\cdots<t_{x}$) the contracted unitary layers collapse into
contractions of only a short portion of the full layer this follows from the
brick property Eq. 65. We demonstrate this process by collapsing a semi-
infinite domain of ‘$+$’ contractions below,

$\raisebox{-0.45pt}{\includegraphics[height=91.04872pt]{Diagrams/contraction-
collapse.png}}.$ (155)

We collapse the ‘$-$’ domain in the same way. Ultimately, every layer reduces
to one of the motifs below,

$\raisebox{-0.45pt}{\includegraphics[height=34.14322pt]{Diagrams/plusminuscontraction.png}}=1+g,\quad\raisebox{-0.45pt}{\includegraphics[height=45.52458pt]{Diagrams/pluszerominus.png}}=\raisebox{-0.45pt}{\includegraphics[height=45.52458pt]{Diagrams/plusperpminus.png}}=-g,\quad\raisebox{-0.45pt}{\includegraphics[height=51.21504pt]{Diagrams/plusperpzerominuscontraction.png}}=g,\quad
M(t)\equiv\raisebox{-0.55pt}{\includegraphics[height=51.21504pt]{Diagrams/Mdef.png}},$
(156)

where $M(t)$ is a single site operator. Each diagram as a whole decomposes
into a products of the motifs below, which we use to introduce a compact
diagrammatic notation. We also give the numerical value of each motif.

$\displaystyle\
\raisebox{-0.45pt}{\includegraphics[height=39.83368pt]{Diagrams/dot-
motif.png}}=-g,\quad\raisebox{-0.2pt}{\includegraphics[height=29.87547pt]{Diagrams/unfilled-
motif.png}}=1+g,$
$\displaystyle\raisebox{-0.45pt}{\includegraphics[height=42.67912pt]{Diagrams/OTOC-
motif.png}}=-g\lambda(t^{\prime}-t),$
$\displaystyle\raisebox{-0.55pt}{\includegraphics[height=56.9055pt]{Diagrams/OTOC-
motif2.png}}=g\lambda(t^{\prime}-t)=-\raisebox{-0.16pt}{\includegraphics[height=28.45274pt]{Diagrams/OTOC-
motif-pic.png}},$ (157)

where we have defined $\lambda(t-t^{\prime})=q\bra{Z^{-}_{t}}M(t+1)\cdots
M(t^{\prime}-1)\ket{Z^{+}_{t^{\prime}}}$. To give an example of the
correspondence between the contracted unitary layer diagrams and this new
diagrammatic notation, consider the diagram corresponding to the sequence
$(t_{1},\cdots,t_{4})=(3,3,6,7)$ below

$\raisebox{-0.45pt}{\includegraphics[height=105.2751pt]{Diagrams/xgreaterthan2_example2.png}}=\raisebox{-0.45pt}{\includegraphics[height=62.59596pt]{Diagrams/xgreaterthan2_example2_shorthand.png}},$
(158)

For $x>2$, the diagrams come in two qualitatively different types: (1)
$t_{1}<t_{2}$ and (2) $t_{1}=t_{2}$. All type 1 diagrams have the following
diagrammatic form

$\raisebox{-0.45pt}{\includegraphics[height=62.59596pt]{Diagrams/type-1_diagram.png}},$
(159)

where $T-t_{x}=t_{1}$. Whereas, type 2 diagrams have the form

$\raisebox{-0.45pt}{\includegraphics[height=62.59596pt]{Diagrams/type-2_diagram.png}}.$
(160)

For $x=2$, we will find diagrams with similar motifs. We had found that for
$x\geq 2$, the relevant diagrams are labelled by a sequence
$(t_{1},t_{2},\cdots,t_{x})$, with $t_{x}=T-t_{1}$ and where $t_{1}=\tau$ is
the length of the OTOCs. For $x=2$ the diagrams are simply labelled by
$(\tau,T-\tau)$. A complication for $x=2$ is the fact that the OTOCs may
overlap in time. This is because no matter where we position the ‘$\perp$’
contraction of site $1$ and ‘$0$’ contraction of site $2$, we can never
encounter any of the vanishing motifs of Eq. D.3, which in the case of $x>2$
force $t_{1}\leq t_{2}$.

To address this complication, we split $x=2$ in three types of diagram: (1)
$T-\tau>\tau$, the OTOCs do not overlap – the treatment of these diagrams is
exactly the same the type-1 diagrams discussed for $x>2$; (2) $\tau=T-\tau$,
the OTOC’s ‘touch’ – this is similar to the type-2 diagrams in $x>2$; (3)
$\tau>T-\tau$, the OTOCs overlap.

The touching OTOC contributions ($\tau=T-\tau$) have the following form/motif,

$\raisebox{-0.45pt}{\includegraphics[height=54.06006pt]{Diagrams/two-OTOC-
motif.png}}=g\lambda(t^{\prime\prime}-t^{\prime})\lambda(t^{\prime}-t).$ (161)

We give the overlapping OTOC contributions ($\tau>T-\tau$) the following
compact notation,

$\raisebox{-0.15pt}{\includegraphics[height=36.98866pt]{Diagrams/overlapping_OTOCs_motif.png}}\
.$ (162)

We study these contributions in detail in the next section, D.4.

We are interested in computing $\int dV\sum_{x}\mathcal{D}^{(2,1)}(x,T)$, we
have seen earlier in this appendix that the $x=0,1$ contributions are
$\order{1/q^{3}}$ (as are the $x<0$ contributions, Eq. 80). Therefore, at
$\order{1/q^{2}}$, we need only sum over $x\geq 2$. Fortunately, there is an
abundance of cancellation between these diagrams. We will cover some examples
and then give the general result. Starting with the simplest, we compute $\int
dV\sum_{x}\mathcal{D}^{(2,1)}(x,T)$ for $T=2,3$ and $4$ explicitly (where the
Haar average is implied, but not written below). The additional factor of
$g^{2}$ in Eq. 144 has been divided through in the equations below.

$\displaystyle\sum_{x}\mathcal{D}^{(2,1)}(x,T=2)/g^{2}$
$\displaystyle=\raisebox{-0.3pt}{\includegraphics[height=31.2982pt]{Diagrams/Dnis1.png}}$
$\displaystyle\sum_{x}\mathcal{D}^{(2,1)}(x,T=3)/g^{2}$
$\displaystyle=\raisebox{-0.38pt}{\includegraphics[height=39.83368pt]{Diagrams/Dnis2.png}}$
$\displaystyle=\raisebox{-0.38pt}{\includegraphics[height=28.45274pt]{Diagrams/Dnis2simplified.png}}$
$\displaystyle\sum_{x}\mathcal{D}^{(2,1)}(x,T=4)/g^{2}$
$\displaystyle=\raisebox{-0.7pt}{\includegraphics[height=79.6678pt]{Diagrams/Dnis3.png}}$
$\displaystyle=\raisebox{-0.28pt}{\includegraphics[height=28.45274pt]{Diagrams/Dnis3simplified.png}}.$

We have used the third rule of Eq. D.3 to cancel the two terms (associated
with $x=2$ and $x=3$) for $T=2$ and to cancel the second and fourth terms and
the fifth and sixth terms for $T=4$. Notice that the terms that remain after
cancellation are all connected diagrams (i.e., the OTOCs either touch or
overlap), all the diagrams with ‘gaps’ (i.e., where a vertical line can be
drawn through them without intersecting a wobbly line, representing an OTOC)
have conspired to cancel. This is no coincidence, it is a consequence of the
fact that processes that contribute to $\Sigma$ (and hence the corrections to
$v_{B}$) must explore only the fast space, this is due to the $Q$ projectors
that project out all slow components at every time-step in $\Sigma$. Diagrams
with a gap represent processes that take a detour to the slow space. We can
see this by returning to the contracted unitary layer picture; take, for
example, the $T=3$ diagram associated with $(t_{1}=1,t_{2}=2)$, this is the
first diagram in the $T=3$ sum above. It is equivalently given by

(163)

Between the two unitary layers we have vectors that are clearly in
$\mathcal{P}$. The only diagrams that have no gaps are those for $x=1$ (which
we have seen all vanish at $\order{1/q^{2}}$) and the overlapping or touching
OTOC diagrams for $x=2$. We will next evaluate the overlapping OTOC diagram
contributions, before finally calculating the touching OTOC diagram
contributions.

###  D.4 Overlapping OTOC diagrams.

In this section we investigate the contributions to
$\mathcal{D}^{(2,1)}(x=2,T)$ that take the form of overlapping OTOCs. We name
this contribution $\mathcal{D}^{(2,1)}_{\textrm{O}}(x=2,T)$. These OTOC
overlap diagrams are given in detail below for OTOC length $\tau$ and total
diagram length $T$.

$\int dV\mathcal{D}^{(2,1)}_{\textrm{O}}(x=2,T)\approx g^{2}\int
dV\sum_{\tau=1}^{T-1}\
\raisebox{-0.55pt}{\includegraphics[height=56.9055pt]{Diagrams/Transfer_matrix_contraction_3.png}}$
(164)

It will be useful to define the single site operators $M_{+}(\tau)$ and
$M_{-}(\tau^{\prime})$ and repeat the definition for $M(t)$ seen in Eq. 156.

$M_{-}(\tau)\equiv\raisebox{-0.55pt}{\includegraphics[height=68.28644pt]{Diagrams/Transfer_matrix_contraction_2.png}},\quad
M_{+}(\tau^{\prime})\equiv\raisebox{-0.55pt}{\includegraphics[height=68.28644pt]{Diagrams/Mplusdef.png}},\quad
M(t)\equiv\raisebox{-0.55pt}{\includegraphics[height=56.9055pt]{Diagrams/Mdef.png}}.$
(165)

The Floquet layers $T-\tau+1$ through to $\tau-1$ have been reduced to two
site operators which has already been introduced in Eq. 103 and named
$\mathcal{T}(t)$ for Floquet layer $t$. With these definitions we have

$\int dV\mathcal{D}^{(2,1)}_{\textrm{O}}(2,T)\approx g^{2}\int
dV\sum_{\tau=1}^{T-1}\
\raisebox{-0.55pt}{\includegraphics[height=56.9055pt]{Diagrams/Transfer_matrix_contraction_4.png}},$
(166)

where we have denoted the following,

.

Using theorem 2 of Sec. 5, the Haar average of a product of two physical OTOCs
is found to $\order{1/q^{2}}$ to be

$\displaystyle\int dV\begin{matrix}\langle
Z\Gamma^{1}_{1}Z(T)\Gamma^{1\dagger}_{\overline{2}}Z\Gamma^{1}_{2}Z(T)\Gamma^{1\dagger}_{\overline{1}}\rangle\\\
\langle
Z\Gamma^{2}_{1}Z(T)\Gamma^{2\dagger}_{\overline{2}}Z\Gamma^{2}_{2}Z(T)\Gamma^{2\dagger}_{\overline{1}}\rangle\end{matrix}=\frac{1}{q^{2}}$
$\displaystyle\left(\delta^{\Gamma^{1}_{1},\Gamma^{2}_{\overline{1}}}\delta^{\Gamma^{1}_{\overline{2}},\Gamma^{2}_{2}}\delta^{\Gamma^{1}_{2},\Gamma^{2}_{\overline{2}}}\delta^{\Gamma^{1}_{\overline{1}},\Gamma^{2}_{1}}\right.$
$\displaystyle\quad+\left.\delta^{\Gamma^{1}_{1},\Gamma^{2}_{\overline{2}}}\delta^{\Gamma^{1}_{\overline{2}},\Gamma^{2}_{1}}\delta^{\Gamma^{1}_{2},\Gamma^{2}_{\overline{1}}}\delta^{\Gamma^{1}_{\overline{1}},\Gamma^{2}_{2}}\right)+\order{1/q^{3}}.$
(167)

Each of the decoration delta constraints can be implemented as described in C.
In doing so, we sandwich each decoration layer with a wiring configuration
labelled $A$ for the first term in Eq. 167 and by a configuration labelled $B$
for the second term. We then use the shorthand below.

$\int
dV\raisebox{-0.55pt}{\includegraphics[height=51.21504pt]{Diagrams/OTOC_product.png}}=\frac{1}{q^{2}}\sum_{a=A,B}\raisebox{-0.55pt}{\includegraphics[height=48.36958pt]{Diagrams/OTOC_product_Haar_av.png}},$
(168)

where the grey boxes are placeholders for the possible decorations at each
layer. A grey box (labeled $t$) may decorated each of the incoming legs $1$
and $2$ with $Z(t)$ and each of the legs $\overline{1}$ and $\overline{2}$
with $Z(t)^{*}$. On the right-hand side, every super-leg carries a label $a$
which labels one of the permutations $A,B$ of the legs
$1,\overline{1},2,\overline{2}$ given in Eq. 169 below,

$\centering\includegraphics[height=25.6073pt]{Diagrams/superleg.png}.\@add@centering$
(169)

Each super-leg is implicitly carrying a factor $1/q^{2}$. The contraction
between decorations within a column is given more explicitly below

$\centering\raisebox{-0.42pt}{\includegraphics[height=34.14322pt]{Diagrams/permutedlegscontraction.png}}=\Tr(R^{a}\Gamma^{\dagger}R^{a}\Gamma^{\prime})/q^{4}=\bra{R^{a}\Gamma
R^{a}}\ket{\Gamma^{\prime}}.\@add@centering$ (170)

All of the onsite scrambling evolution can be dropped on the right-hand-side
of Eq. 168 as the leg contractions are between decorations at the same time.
In the case of Eq. 166, the OTOCs are off-set by $\tau^{\prime}=T-\tau$
Floquet time-steps. These OTOCs can be brought into alignment by globally
shift the time arguments in the OTOCs, then we can use Eq. 168. After Haar
averaging, and shifting the OTOCs back to their original positions, the leg
contractions will stretch over $\tau^{\prime}$ steps, as shown below.

$\frac{g^{2}}{q^{2}}\sum_{s\geq 0}\sum_{t\geq
0}\delta^{s+2t+3,T}\sum_{a=A,B}\raisebox{-0.55pt}{\includegraphics[height=56.9055pt]{Diagrams/Transfer_matrix_contraction_5.png}},$
(171)

where all scrambling evolution time arguments have been dropped. The new
variables are related to those in Eq. 166 by $t=T-\tau-1\geq 0$ and
$s=2\tau-T-1\geq 0$. As an example, we show the $t=1$, $s=2$ case,

$\raisebox{-0.4pt}{\includegraphics[height=39.83368pt]{Diagrams/Transfer_matrix_contraction_5_ex.png}}=\raisebox{-0.4pt}{\includegraphics[height=36.98866pt]{Diagrams/Transfer_matrix_chain_ex1.png}},$
(172)

where on the right hand side we have introduced a shorthand which makes
obvious the decomposition into two chains of tensor contractions. In general,
the Haar average decomposes into $t+1$ chains. By summing over $T\geq 1$, we
are able to drop the delta constraint above, doing this sum is equivalent to
calculating the Laplace transformed $\int dV\mathcal{D}^{2,1}(x=2,T)$ at
$z=0$, i,e $\int dV\mathcal{D}^{2,1}(x=2,z=0)$. This is sufficient for
calculating $v_{B}$. Before we determine the decomposition for a general
$s\geq 0$ and $t\geq 0$, we first define the following chains,

$\displaystyle\raisebox{-0.4pt}{\includegraphics[height=45.52458pt]{Diagrams/Chain.png}},\quad\raisebox{-0.4pt}{\includegraphics[height=45.52458pt]{Diagrams/Chain+.png}},$
$\displaystyle\raisebox{-0.4pt}{\includegraphics[height=45.52458pt]{Diagrams/Chain-.png}},\quad\raisebox{-0.4pt}{\includegraphics[height=45.52458pt]{Diagrams/Chain+-.png}}.$
(173)

For any $(s,t)$ There are two qualitatively distinct types of contribution to
Eq. 171: (1) $t+1\equiv 0\pmod{s+1}$, a product of $t$ $C^{a}$ chains and one
$C^{a}_{+-}$ chain; (2) otherwise, a product of $t-1$ $C^{a}$ chains, one
$C^{a}_{+}$ chain and one $C^{a}_{-}$ chain. For case 1, the $C^{a}_{+-}$
chain is $n=\frac{s+1}{t+1}-1\geq 0$ transfer matrices long, i.e.,
$C^{a}_{+-}(n)$, while all $t$ of $C^{a}$ chains are $n+1$ matrices long,
$C^{a}(n+1)$. For case 2 with $s+1=n(t+1)+k$ for $1\leq k\leq t$ and starting
from the left, the first $k-1$ $M$’s are on chains of length $n+1$ and
terminate on an $M$. The $k$-th $M$ sits on a chain of length $n$ and
terminates on $M_{-}$, the following $t-k$ $M$’s sit on chains length $n$ that
terminate on an $M$. Finally, the $M_{+}$ sits on a chain of length $n$ and
terminates on an $M$. All together, the contribution is
$C^{a}_{+}(n)C^{a}_{-}(n)C^{a}(n+1)^{k-1}C^{a}(n)^{t-k}$. This is summarised
below,

$\frac{g^{2}}{q^{2}}\sum_{a=A,B}\sum_{n\geq 0}\left[\sum_{t\geq
0}C^{a}_{+-}(n)C^{a}(n+1)^{t}+\sum_{t\geq
1}\sum_{k=1}^{t}C^{a}_{+}(n)C^{a}_{-}(n)C^{a}(n+1)^{k-1}C^{a}(n)^{t-k}\right].$
(174)

In both cases, all but the $n$ sums can be evaluated to give,

$\displaystyle\int dV\mathcal{D}^{(2,1)}_{\textrm{O}}$
$\displaystyle(x=2,z=0)=$
$\displaystyle\frac{g^{2}}{q^{2}}\sum_{a=A,B}\sum_{n\geq
0}\left[\frac{C^{a}_{+-}(n)}{1-C^{a}(n+1)}+\frac{C^{a}_{+}(n)C^{a}_{-}(n)}{(1-C^{a}(n))(1-C^{a}(n+1))}\right].$
(175)

It remains to calculate the different chains. To do this we approach the
problem as a transfer matrix problem, where $\mathcal{T}$ is the transfer
matrix. Using the definition of $\mathcal{T}$ in Eq. 103, the decoration
decomposition of a two-site brick in Eq. 62 and definitions of $M$ and
$M^{\pm}$ in Eq. 165, $\mathcal{T}$ and each of $M$, $M_{+}$ and $M_{-}$ have
the following properties,

$\raisebox{-0.4pt}{\includegraphics[height=29.87547pt]{Diagrams/transfer_matrix_properties.png}},\quad\raisebox{-0.4pt}{\includegraphics[height=17.07182pt]{Diagrams/M_properties.png}}.$
(176)

Then, using the definition of the labelled super legs in Eq. 169, the chain
$C^{a}(n)$ simplifies to the expression below,

$C^{a}(n)=\raisebox{-0.4pt}{\includegraphics[height=36.98866pt]{Diagrams/Transfer_matrix_simple_chain.png}}.$
(177)

Notice that the dependence on the label $a$ vanished. This is true of all
chains. Algebraically, these chains (now without leg labels) are equivalently
given by

$\displaystyle C(n)=\bra{M}\mathcal{T}^{n}\ket{M}$ $\displaystyle,\quad
C_{+}(n)=\bra{M_{+}}\mathcal{T}^{n}\ket{M},$ $\displaystyle
C_{+}(n)=\bra{M}\mathcal{T}^{n}\ket{M_{-}}$ $\displaystyle,\quad
C_{+-}(n)=\bra{M_{+}}\mathcal{T}^{n}\ket{M_{-}},$ (178)

where these angles braces reflect the trace inner product for tensors with
input and output super legs $l=(1,\overline{1},2,\overline{2})$.

$\bra{B}\ket{A}=\Tr(B^{\dagger}A)/q^{4}=\raisebox{-0.4pt}{\includegraphics[height=28.45274pt]{Diagrams/Trace_inner_prod_super_leg.png}}.$
(179)

Equipped with this, we can now write, for the contribution due to overlapping
OTOC diagrams, the following

$\int
dV\mathcal{D}^{(2,1)}_{\textrm{O}}(x=2,z=0)=\frac{2g^{2}}{q^{2}}f(\varepsilon),$
(180)

where $f(\varepsilon)$ is given by

$f(\varepsilon)=\sum_{n\geq
0}\left[\frac{C_{+-}(n)}{1-C(n+1)}+\frac{C_{+}(n)C_{-}(n)}{(1-C(n))(1-C(n+1))}\right],$
(181)

where the chains are implicitly dependent on $\varepsilon$. To evaluate this
sum, we must understand the space that the transfer matrix acts on. Each of
the legs $1$, $\overline{1}$, $2$ and $\overline{2}$ may be either undecorated
or carry a $Z$ decoration. This means that the state our state space is
dimension $2^{4}$ and $\mathcal{T}$ is a $16\times 16$ matrix. Because all the
$M_{(\pm)}$ are even in the number of $Z$ decorations, and $\mathcal{T}$
preserves decoration parity, we are able to reduce the state space to those
states with an even number of $Z$ decorations only, i.e., $8$ states. It will
be useful to use the basis below,

$\raisebox{-0.4pt}{\includegraphics[height=65.44142pt]{Diagrams/Basis.png}}.$
(182)

Let $\mathcal{S}$ be the space spanning
$\\{\ket{\mathbb{1}},\ket{S_{1\overline{1}}},\ket{S_{1\overline{2}}},\ket{S_{12}},\ket{4}\\}$
and $\mathcal{A}$ be the space spanning
$\\{\ket{A_{1\overline{1}}},\ket{A_{1\overline{2}}},\ket{A_{12}}\\}$ One can
easily check that $\ket{M}$, $\ket{M_{+}}$ and $\ket{M_{-}}$ are all in
$\mathcal{S}$. Explicitly, with $u(\varepsilon)=\sin(\varepsilon)^{2}$

$\displaystyle\ket{M}=(u^{4}+(1-u)^{4})\ket{\mathbb{1}}+2u^{2}(1-u)^{2}\ket{4}-(1+g)\frac{g}{\sqrt{2}}\left(\ket{S_{1\overline{1}}}+\ket{S_{1\overline{2}}}\right)+\frac{g^{2}}{\sqrt{2}}\ket{S_{12}},$
(183)
$\displaystyle\ket{M_{+}}=-\frac{g(1+g)}{\sqrt{2}}\left(\ket{S_{12}}-\ket{S_{1\overline{2}}}\right)+\frac{g^{2}}{2}\left(\sqrt{2}\ket{S_{1\overline{2}}}-\ket{\mathbb{1}}-\ket{4}\right),$
(184)
$\displaystyle\ket{M_{-}}=-\frac{g(1+g)}{\sqrt{2}}\left(\ket{S_{12}}-\ket{S_{1\overline{1}}}\right)+\frac{g^{2}}{2}\left(\sqrt{2}\ket{S_{1\overline{1}}}-\ket{\mathbb{1}}-\ket{4}\right).$
(185)

$\mathcal{T}$ is block diagonal in the subspaces $\mathcal{S}$ and
$\mathcal{A}$, and since all $\ket{M_{(\pm)}}$ lie in $\mathcal{S}$, we are
able to restrict our considerations to this five dimensional space only. In
this restricted space with basis order
$(\mathbb{1},4,S_{1,\overline{1}},S_{1,\overline{2}},S_{1,2})$, $\mathcal{T}$
is given by the product $\mathcal{T}=\tilde{T}_{-}\tilde{U}\tilde{T}_{+}$,
where

$\tilde{U}=\left(\begin{array}[]{ccccc}(1-u)^{2}&0&0&0&0\\\ 0&u^{2}&0&0&0\\\
0&0&-\frac{g}{2}&0&0\\\ 0&0&0&-\frac{g}{2}&0\\\ 0&0&0&0&\frac{g}{2}\\\
\end{array}\right)$ (186)
$\tilde{T}_{-}=\left(\begin{array}[]{ccccc}(1-u)^{2}&u^{2}&-\frac{g}{\sqrt{2}}&0&0\\\
u^{2}&(1-u)^{2}&-\frac{g}{\sqrt{2}}&0&0\\\
-\frac{g}{\sqrt{2}}&-\frac{g}{\sqrt{2}}&g+1&0&0\\\ 0&0&0&g+1&-g\\\
0&0&0&-g&g+1\\\ \end{array}\right)$ (187)
$\tilde{T}_{+}=\left(\begin{array}[]{ccccc}(1-u)^{2}&u^{2}&0&-\frac{g}{\sqrt{2}}&0\\\
u^{2}&(1-u)^{2}&0&-\frac{g}{\sqrt{2}}&0\\\ 0&0&g+1&0&-g\\\
-\frac{g}{\sqrt{2}}&-\frac{g}{\sqrt{2}}&0&g+1&0\\\ 0&0&-g&0&g+1\\\
\end{array}\right)$ (188)

Four of the eigenvalues of $\mathcal{T}$ are bounded by
$\absolutevalue{\lambda_{1,2,3,4}}\leq\absolutevalue{g(\varepsilon)}(1-2\absolutevalue{g(\varepsilon)})/2$.
The largest eigenvalue (for all $\varepsilon$) is bounded by
$\absolutevalue{\lambda_{5}}\leq(1-\absolutevalue{g(\varepsilon)})^{3}$.

In Sec. 4 we showed that $\langle v_{B}(\varepsilon)\rangle$ must have the
symmetries $\varepsilon\to-\varepsilon$ and $\varepsilon\to\pi/2+\varepsilon$.
All analytically computed $\order{1/q^{2}}$ contributions (the $(a,b)=(4,4)$
contribution and the touching OTOC contribution) are found to respect this
symmetry, as does the $\order{1}$ contribution from $\Omega$. We therefore
conclude that $f(\varepsilon)$ must also have this symmetry. Additionally, we
know that $f(\varepsilon)$ is a function of $\sin(\varepsilon)^{2}$ only (the
transfer matrix and the initial and final vectors in our transfer matrix
calculation are functions of $\sin(\varepsilon)^{2}$ only). Together with the
symmetry requirements, this means that $f(\varepsilon)$ is in fact a function
of $s(\varepsilon)=\sin(\varepsilon)^{2}\cos(\varepsilon)^{2}$.

We evaluate the sum in Eq. 181 analytically for small $\varepsilon$, finding
$f(\varepsilon)\approx\varepsilon^{2}/7$. Knowing that $f$ is a function of
$s(\varepsilon)$ only, we use this to factorise
$f(\varepsilon)=\frac{1}{7}s(\varepsilon)w(s(\varepsilon))$, where
$w(\varepsilon)$ is a function that approaches $1$ as $\varepsilon\to 0$. The
transfer matrix dramatically simplifies at the point $s(\varepsilon)=1/4$.
Where all but one eigenvalue is zero. We find analytically that
$f(\varepsilon)\propto(1-4s)^{2}$ around this point. This suggests a further
factorisation
$f(\varepsilon)=\frac{1}{7}s(\varepsilon)(1-4s)^{2}p(s(\varepsilon))$. We find
a very good quadratic polynomial approximation for
$p(\varepsilon)=1+as+bs^{2}$, with $a=6.8$ and $b=16.1$.

$f(\varepsilon)=\frac{1}{7}s(\varepsilon)(1-4s(\varepsilon))^{2}(1+as(\varepsilon)+bs(\varepsilon)^{2})$
(189)

.

Figure 5: The contribution to $v_{B}$ from processes with overlapping OTOCs is
given by $2g(\varepsilon)^{2}f(\varepsilon)/q^{2}$, $f(\varepsilon)$ is by
numerically and plotted in blue. A fitting function is plotted in red.

We also find numerically that $\int dV\mathcal{D}^{2,1}_{\textrm{O}}(x=2,T)$
decays exponentially quickly, with a decay rate bounded below by
$\gamma_{\textrm{O}}(\varepsilon)\geq 6s(\varepsilon)$.

###  D.5 Touching OTOC diagrams

We will now sum up all of the touching OTOC diagrams. These are given in Eq.
161 with $T=2\tau$ so that the two OTOC have the same length. We name this
contribution $\mathcal{D}^{(2,1)}_{\textrm{Touch}}(x=2,T)$

We are only interested in the $z\to i0^{+}$ limit, this is given by

$\int dV\mathcal{D}^{(2,1)}_{\textrm{Touch}}(k=0,z=0)=g^{2}\int
dV\sum_{t=1}^{\infty}g\lambda(t)^{2}+\order{1/q^{3}}.$ (190)

The factor $g\lambda(t)^{2}$ is the contribution from the diagrams Eq. 161.
The Haar average of a product of OTOCs given in Eq. 168 is deployed again to
say that

$\int dV\lambda(t)^{2}=\frac{2}{q^{2}}\bra{M}\ket{M}^{t-1}.$ (191)

This gives

$\int
dV\sum_{t=1}^{\infty}g^{3}\lambda(t)^{2}=\frac{2g^{3}}{q^{2}}\frac{1}{1-\bra{M}\ket{M}}=-\frac{2g^{2}}{q^{2}}\nu(\varepsilon),$
(192)

where $\bra{M}\ket{M}=1-8s(1-2s)(1-s(1-2s))$,
$s(\varepsilon)=\sin(\varepsilon)^{2}\cos(\varepsilon)^{2}$ and
$\nu(\varepsilon)=-g/(1-\bra{M}\ket{M})$. Touching OTOC diagrams decay
exponentially with a decay rate
$\gamma_{\textrm{Touch}}(\varepsilon)=\log(\bra{M}\ket{M}^{-1})\geq
8s(\varepsilon)$.

###  D.6 Summary

Recalling Eq. 80, the contributions from $x<0$ are $\order{1/q^{4}}$ or
smaller. Piecing together the $\order{1/q^{2}}$ corrections, we find $\int
dV\mathcal{D}^{(2,1)}(k=0,z=0)$ is given by

$\int
dV\mathcal{D}^{(2,1)}(k=0,z=0)=\frac{2g^{2}}{q^{2}}(f(\varepsilon)-\nu(\varepsilon))+\order{1/q^{3}}.$
(193)

##  Appendix E Remaining $(a,b)$

###  E.1 $(1,2)$

$\mathcal{D}^{1,2}_{\Gamma}(x\geq 0,T)=g^{2}\times\left[\
\begin{matrix}\vspace{-0.15cm}\\\
{\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>}0}\\\
{\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>}1}\\\ \\\ \\\ \\\ \\\ \end{matrix}\ \ \begin{matrix}\ \ \vdots\\\ \ \
\bra{\phi_{+}}\\\ q\bra{-}\\\ \ \ \vdots\\\ q\bra{-}\\\ \ \ \bra{-}\\\ \ \
\vdots\end{matrix}\quad\framebox{ $\begin{matrix}\\\ \\\ \\\ \Gamma\\\ \\\ \\\
\\\ \end{matrix}$}\quad\begin{matrix}\vdots\\\ \ket{+}\\\ \ket{+}\\\ \vdots\\\
\ket{+}\\\ \ket{\phi_{-}(T)}\\\ \vdots\end{matrix}\ \ \begin{matrix}\\\ \\\
\\\ \\\ \vspace{-1.01cm}\\\
{\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>}x}\\\
{\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>}x+1}\end{matrix}\right]$ (194)

Here, we can simply use Eq. 89 to see that both sites $0$ and $x+1$ contribute
a product of non-trivial correlation functions plus terms of size $1/q^{2}$.
The Haar average of this is $\order{1/q^{4}}$.

$\int dV\mathcal{D}^{1,2}_{\Gamma}(x\geq 0,T)=\order{1/q^{4}}.$ (195)

###  E.2 $(3,3)$

  * •

$x=0$:

$\mathcal{D}^{3,3}_{\Gamma}(x=0,T)=g^{2}\times\left[\
\begin{matrix}{\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>}0}\\\
{\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>}1}\end{matrix}\ \ \begin{matrix}\bra{+}Z^{\otimes 2}\\\ \bra{-}Z^{\otimes
2}\end{matrix}\ \Gamma\ \begin{matrix}Z(T)^{\otimes 2}\ket{+}\\\ Z(T)^{\otimes
2}\ket{-}\end{matrix}\right]\left(\prod_{r<0}\bra{+}\Gamma^{r}\ket{+}\right)\left(\prod_{r>1}\bra{-}\Gamma^{r}\ket{-}\right)$
(196)

In this case, both contours on site $0$ and on site $1$ are non-trivially
decorated, contributing a total of four non-trivial correlators. The Haar
average is then $\order{1/q^{4}}$.

  * •

$x\geq 1$: In this case, each site $0$ and $x+1$ contribute a product of non-
trivial correlators and therefore, as above, are $\order{1/q^{4}}$ or smaller.

All together

$\int dV\mathcal{D}^{1,2}_{\Gamma}(x\geq 0,T)=\order{1/q^{4}}.$ (197)

###  E.3 $(3,4)$ and $(4,3)$

  * •

$x=0$:

$\mathcal{D}^{3,4}_{\Gamma}(0,T)=-ihg\times\left[\
\begin{matrix}{\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>}0}\\\
{\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>}1}\end{matrix}\ \ \begin{matrix}\bra{+}Z^{\otimes 2}\Gamma^{0}K(T)\ket{+}\\\
\bra{-}Z^{\otimes
2}\Gamma^{1}K(T)\ket{-}\end{matrix}\right]\left(\prod_{r<0}\bra{+}\Gamma^{r}\ket{+}\right)\left(\prod_{r>1}\bra{-}\Gamma^{r}\ket{-}\right)$
(198)

Both site $0$ and site $1$ contribute two non-trivial correlators. The Haar
average of this term is $\order{1/q^{4}}$ or smaller.

  * •

$x\geq 1$: Each of the sites $0$ and $x+1$ contribute a product of two non-
trivial correlators. Each of the sites $1$ and $x$ contribute a non-trivial
correlator. If $x=1$, then this is only five non-trivial correlators,
otherwise it is six. The Haar average is then $\order{1/q^{5}}$

All together,

$\int dV\mathcal{D}^{4,4}_{\Gamma}(x\geq 0,T)=\order{1/q^{4}}.$ (199)

###  E.4 $(a,b)$, where either $a\in{1,2}$ and $b\in{3,4}$ or the converse

In these case, Eq. 89 is again enough to show that $\int
dV\mathcal{D}^{a,b}_{\Gamma}(x,T)=\order{1/q^{3}}$.

##  Appendix F Model variation: Independent scramblers $V_{x,t}$

We can carry out a similar calculation for a version of the model with
independently distributed scramblers $V_{x,t}$ for each site and at each layer
of unitaries $U_{t}$ (analogous to the Floquet layer in the Floquet model),
i.e., we break both spatial and temporal translation symmetry. We start with,
for large $t$, the fact $\sum_{x}x\rho(x,t)=v_{B}t$ to write

$\displaystyle v_{B}$
$\displaystyle=\lim_{t\to\infty}\sum_{x}x\left[\rho(x,t)-\rho(x,t-1)\right]=\lim_{t\to\infty}\lim_{k\to
0}\sum_{x}\mathrm{i}\partial_{k}e^{-\mathrm{i}kx}\left[\rho(x,t)-\rho(x,t-1)\right]$
$\displaystyle=\lim_{t\to\infty}\lim_{k\to
0}\mathrm{i}\partial_{k}\left[\rho(k,t)-\rho(k,t-1)\right].$ (200)

Using the definition of $\rho(k,n)$ we write the following

$\displaystyle\rho(k,n)$ $\displaystyle=\left(W^{k}|U_{0}(P+Q)U_{1}\cdots
U_{n}|W^{k}\right)=\left(W^{k}|U_{0}|W^{k}\right)\rho(k,n-1)+H(k,n)$
$\displaystyle=(1+\Omega(k))\rho(k,n-1)+H(k,n)$ (201)

where $H(k,n)$ is given by

$\displaystyle H(k,n)$ $\displaystyle=\left(W^{k}|U_{0}QU_{1}\cdots
U_{n}|W^{k}\right)$ $\displaystyle=\left(W^{k}|U_{0}QU_{1}\cdots
U_{n-1}(P+Q)U_{n}|W^{k}\right)$
$\displaystyle=(1+\Omega(k))H(k,n-1)+\sigma(k,n),$ (202)

where $\sigma(k,n)$ is as previously defined, but with the independently
distributed scramblers.

$\sigma(k,n)=\left(W^{k}|L_{0}QU_{1}\cdots U_{n-1}QL_{n}|W^{k}\right).$ (203)

Using Eq. F inductively and the fact that $H(k,n=0)=0$ we find

$\rho(k,n)-\rho(k,n-1)=\Omega(k)\rho(k,n-1)+\sum_{m=1}^{n}(1+\Omega(k))^{n-m-1}\sigma(k,m).$
(204)

Using this in Eq. F, we find

$v_{B}=\lim_{n\to\infty}\lim_{k\to
0}\mathrm{i}\partial_{k}\left[\Omega(k)\rho(k,n-1)+\sum_{m=1}^{n}(1+\Omega(k))^{n-m-1}\sigma(k,m)\right].$
(205)

Then writing $\Omega(k)=v_{0}(1-e^{-\mathrm{i}k})$ where
$v_{0}=\absolutevalue{g(\varepsilon)}$ is the $q\to\infty$ butterfly velocity,
and noticing that $\lim_{k\to
0}\partial_{k}\left(\Omega(k)\rho(k,n)\right)=\mathrm{i}v_{0}$, we find

$\delta v_{B}\equiv v_{B}-v_{0}=\lim_{n\to\infty}\lim_{k\to
0}\mathrm{i}\partial_{k}\sum_{m=1}^{n}(1+\Omega(k))^{n-m-1}\sigma(k,m).$ (206)

Analogous to the definition in Eq. 76 for the Floquet model, we define the
quantity,

$\mathcal{D}(x\to y,n)=q^{y-x}\bra{F^{x}}L_{0}QU_{1}\cdots
U_{n-1}QL_{n}\ket{F^{y}}$ (207)

then $\sigma(k,n)$ is given by

$\displaystyle\sigma(k,n)=\frac{1}{L(1-q^{-2})}\sum_{x,y}$ $\displaystyle
e^{\mathrm{i}k(x-y)}\bigg{[}\mathcal{D}(x\to y,n)-\mathcal{D}(x\to y-1,n)$
$\displaystyle\left.-\frac{1}{q^{2}}\left(\mathcal{D}(x-1\to
y,n)-\mathcal{D}(x-1\to y-1,n)\right)\right]$ (208)

The circuit averaged $\mathcal{D}(x\to y,n)$ is translationally invariant
$\langle\mathcal{D}(x\to y,n)\rangle=\langle\mathcal{D}(y-x,n)\rangle$ where
$\mathcal{D}(x,n)=\mathcal{D}(0\to x,n)$. The circuit average of $\sigma(k,n)$
is given by

$\langle\sigma(k,n)\rangle=\eta(k)(1-e^{-\mathrm{i}k})\sum_{x}e^{-\mathrm{i}kx}\langle\mathcal{D}(x,n)\rangle=\eta(k)(1-e^{-\mathrm{i}k})\langle\mathcal{D}(k,n)\rangle.$
(209)

The circuit averaging of Eq. 206 is given by

$\displaystyle\langle\delta v_{B}\rangle$
$\displaystyle=\lim_{n\to\infty}\lim_{k\to
0}\mathrm{i}\partial_{k}\sum_{m=1}^{n}(1+\Omega(k))^{n-m-1}\eta(k)(1-e^{-\mathrm{i}k})\langle\mathcal{D}(k,m)\rangle$
$\displaystyle=-\sum_{m=1}^{\infty}\langle\mathcal{D}(k=0,m)\rangle=-\langle\mathcal{D}(k=0,z=0)\rangle,$
(210)

This is same expression for the circuit averaged correction to $v_{B}$ that we
found using the memory matrix formalism in Eq. 112. We calculate
$\langle\mathcal{D}(k=0,z=0)\rangle$ in the same way as in Sec. 6. However,
the two $\order{1/q^{2}}$ contributions found for the Floquet model (i.e., the
so-called $(2,1)$ and $(4,4)$ terms) relied on correlations between different
sites and on the discrete time translation symmetry. With independently
distributed scramblers, such terms are smaller than $\order{1/q^{2}}$. This
gives the result $\langle\delta v_{B}\rangle=\order{1/q^{3}}$ which leads to

$\langle v_{B}\rangle=\frac{1-\cos(4\varepsilon)}{4}+\order{1/q^{3}}.$ (211)

We can employ the same calculation above for the variant of the model with
spatial translation symmetry but no time translation symmetry (i.e.,
independently random scramblers between time-steps) in order to find
$\langle\delta v_{B}\rangle=\delta v_{S}(\varepsilon)+\order{1/q^{3}}$ where
$\delta v_{S}(\varepsilon)$ is as given in Eq. 113.

$\langle v_{B}\rangle=\frac{1-\cos(4\varepsilon)}{4}+\delta
v_{S}(\varepsilon)+\order{1/q^{3}},\quad\delta
v_{S}(\varepsilon)=\frac{1}{q^{2}}\frac{1+5s-4s^{2}}{1-s-3s^{2}}.$ (212)

Comparing these result to the result of Eq. 112 allows us to see that the
$\order{1/q^{2}}$ corrections only arise when there is spatial translation
symmetry. Moreover, this enables us to identify the $\delta v_{F}$
contribution (see Eq. 113) with the presence of time translation symmetry.

###  F.1 The $\varepsilon\to 0$ limit in the absence of time-translation
symmetry

In Sec. 6.2.2, we identified all the contributions to $\sigma(x=0,t)$ that
contained the product of two non-trivial correlators and used theorem 2 to
evaluate these contributions to $\order{1/q^{2}}$. We were unable to extract
the correct $\varepsilon\to 0$ behaviour of these contributions due to the
appearance of $\varepsilon$ in the denominator (after summing over the time
$t$). Luckily, for the model variant with independently random scramblers at
each time-step, we are able to accurately sum over all ‘two-correlator’
contributions, after circuit averaging. While this is ignores correlations
between Floquet layers, it demonstrated the failure of naively expanding the
memory matrix in $1/q$ before taking the $\varepsilon\to 0$ limit and the
resolution by summing over all contributions.

We identify the contributions that have the form of two non-trivial
correlators in the same way as done in Sec. 6.2.2. We repeat this below for
ease of reading.

$\mathcal{D}^{4,4}(x=0,T)=-h^{2}\sum_{i,j\in\\{1,2\\}}\begin{matrix}\bra{+}Z_{i}\\\
\bra{-}Z_{j}\end{matrix}\ \framebox{ $\begin{matrix}\vspace{-5mm}\\\
1\vspace{1mm}\end{matrix}$}\ \framebox{ $\begin{matrix}\vspace{-5mm}\\\
2\vspace{1mm}\end{matrix}$}\ \cdots\ \framebox{
$\begin{matrix}\vspace{-5mm}\\\ T-1\vspace{1mm}\end{matrix}$}\
\begin{matrix}Z(T)_{i}\ket{+}\\\ Z(T)_{j}\ket{-}\end{matrix}\ +\ \cdots,$
(213)

where, so as to ensure that only sites $0$ and $1$ contribute non-trivial
correlators, we have sandwiched each layer by $\bra{+}$ ($\bra{-}$) and
$\ket{+}$ ($\ket{-}$) for sites $x<0$ ($x>1$) to leave the two site operators
$\mathcal{T}(t)$ introduced in Eq. 103 and represented above in Eq. 213 by the
rectangle blocks. The second ellipsis represents terms that contain at least
three non-trivial correlators. Using the same argument as was used in Sec.
6.2.2, we see that the $(i,j)=(1,1)$ and $(2,2)$ choices give identical
contributions and the $(1,2)$ and $(2,1)$ contributions are the complex
conjugate of the $(i,i)$ contributions, Eq. 105. Therefore, we are free to
study only the $(i,j)=(1,1)$ case. We insert projectors ensure that the
$2,\overline{1}$ contour of site $0$ and the $2,\overline{2}$ contour of site
$1$ remain undecorated (so as to ensure that there are only two non-trivial
correlators). The contribution is given below,

$\raisebox{-0.45pt}{\includegraphics[height=91.04872pt]{Diagrams/44-diagram.png}},$
(214)

where each closed loops is associated with a normalising factor of $1/q$ that
is suppressed in the diagram. The contracted single time-step unitaries are
now just operators on four indices. We denote the contribution to
$\mathcal{D}^{4,4}(x=0,T)$ that includes only two non-trivial contours as
$\mathcal{D}_{2}^{4,4}(x=0,T)$. Separating out the scrambling part of the
unitaries we find

$\mathcal{D}_{2}^{4,4}(x=0,T)=-\frac{2h^{2}}{q^{2}}\
\raisebox{-0.45pt}{\includegraphics[height=56.9055pt]{Diagrams/44_independent_scramblers.png}}+\textrm{c.c},$
(215)

where each of the pink bricks is given by the following contraction of
$\mathcal{T}$.

$Y=\raisebox{-0.45pt}{\includegraphics[height=71.13188pt]{Diagrams/contracted_brick.png}}\times\frac{1}{q^{2}}.$
(216)

Circuit averaging this tensor contraction is easy as each unitary only appears
twice. We use the following result for the Haar average of two $V$’s and two
$V^{\dagger}$’s.

$\int
dV\raisebox{-0.5pt}{\includegraphics[height=56.9055pt]{Diagrams/two_U_Haar_integral.png}}=K,$
(217)

where $K$ is a projector. Using the definitions of $\ket{+}$, $\ket{-}$,
$\ket{0}$ and $\ket{\perp}$, $K$ is given algebraically as

$K\equiv\frac{\ket{-}\bra{\perp}}{1-q^{-2}}+\frac{\ket{+}\bra{0}}{1-q^{-2}}=\frac{\ket{\perp}\bra{-}}{1-q^{-2}}+\frac{\ket{0}\bra{+}}{1-q^{-2}}.$
(218)

Using this, we write

$\displaystyle\langle\mathcal{D}_{2}^{4,4}(x=0,T)\rangle$
$\displaystyle=-2h^{2}\bra{-}Z^{\otimes 2}(KYK)^{T-1}Z^{\otimes
2}\ket{-}+\textrm{c.c}$
$\displaystyle=-\frac{2h^{2}}{q^{2}}\bra{0}(KYK)^{T-1}\ket{0}+\textrm{c.c}.$
(219)

The matrix $KYK$ is a $2\times 2$ matrix acting on the space of wirings
$\textrm{Span}\\{\ket{+},\ket{-}\\}$. Choosing the orthormal basis
$\ket{S,A}=\frac{1}{\sqrt{2}}(\ket{+}\pm\ket{-})$. $KYK$ is given by

$KYK=\frac{\xi}{2}\left(\begin{array}[]{cc}1&1\\\
1&1\end{array}\right)+\frac{1+g}{2}\left(\begin{array}[]{cc}1&-1\\\
-1&1\end{array}\right)+\frac{1+g+ih}{q}\left(\begin{array}[]{cc}1&0\\\
0&-1\end{array}\right).$ (220)

By diagonalising this, we find

$\sum_{n=0}^{\infty}\frac{1}{q^{2}}\bra{0}(KYK)^{n}\ket{0}=\frac{2}{1+8q^{2}\varepsilon^{2}+\order{q^{2}\varepsilon^{3},q\varepsilon^{2}}}.$
(221)

Using $h(\varepsilon)\approx\varepsilon$ for small $\varepsilon$,
$\langle\mathcal{D}_{2}^{4,4}(x=0,z=0)\rangle$ is then given by

$\langle\mathcal{D}_{2}^{4,4}(x=0,z=0)\rangle=-\frac{8\varepsilon^{2}}{1+8q^{2}\varepsilon^{2}+\order{q^{2}\varepsilon^{3},q\varepsilon^{2}}}.$
(222)

Importantly, this contribution, which keeps account of all orders in $q$,
vanishes as $\varepsilon\to 0$, unlike the $\order{1/q^{2}}$ contribution
(naively) calculated in Sec. 6.2.2, which (incorrectly) predicts
$\lim_{\varepsilon\to 0}\langle\mathcal{D}^{4,4}(x=0,z=0)\rangle=-1/q^{2}$.
However, this is still far from a complete account of contributions to
$\langle\mathcal{D}(k=0,z=0)\rangle$ (for instance, contributions involving
more than two correlators), but it does demonstrate that the apparent
pathological behaviour as $\varepsilon\to 0$ is the result of a naive (and
incorrect at small $\varepsilon$) counting of the $\order{1/q^{2}}$
contributions. Notice that for $q\varepsilon\gg 1$, the small $\varepsilon$
results of Sec. 6.2.2 (which predict
$\langle\mathcal{D}^{4,4}(x=0,z=0)\rangle\approx-1/q^{2}$) are in agreement
with the analysis presented here.