# Supplementary Material to “Evidence for bootstrap percolation dynamics in a
photoinduced phase transition”

Tyler Carbin<EMAIL_ADDRESS>Department of Physics and Astronomy, UCLA,
Los Angeles, CA 90095-1547  Xinshu Zhang  Department of Physics and Astronomy,
UCLA, Los Angeles, CA 90095-1547  Adrian B. Culver  Department of Physics and
Astronomy, UCLA, Los Angeles, CA 90095-1547  Mani L. Bhaumik Institute for
Theoretical Physics, Department of Physics and Astronomy, University of
California Los Angeles, Los Angeles, CA 90095  Hengdi Zhao  Department of
Physics, University of Colorado at Boulder, Boulder, Colorado 80309, USA
Alfred Zong  Department of Chemistry, University of California at Berkeley,
Berkeley, CA, 94720, USA  Materials Sciences Division, Lawrence Berkeley
National Laboratory, Berkeley, CA, 94720, USA  Rishi Acharya  Department of
Physics and Astronomy, UCLA, Los Angeles, CA 90095-1547  Cecilia J. Abbamonte
Department of Physics and Astronomy, UCLA, Los Angeles, CA 90095-1547  Rahul
Roy  Department of Physics and Astronomy, UCLA, Los Angeles, CA 90095-1547
Mani L. Bhaumik Institute for Theoretical Physics, Department of Physics and
Astronomy, University of California Los Angeles, Los Angeles, CA 90095  Gang
Cao  Department of Physics, University of Colorado at Boulder, Boulder,
Colorado 80309, USA  Anshul Kogar  Department of Physics and Astronomy, UCLA,
Los Angeles, CA 90095-1547<EMAIL_ADDRESS>
###### Contents

  1. I Rotational-Anisotropy SHG fitting
    1. I.1 Imperfect alignment
    2. I.2 Multiple Domains
  2. II SHG Wavelength Dependence
  3. III Time Trace Fits
  4. IV Recovery time traces
  5. V Estimates of Time Resolution
  6. VI Calculation of $d_{\gamma}$
  7. VII Penetration depth calculations
  8. VIII Considerations of homogeneous mechanisms
  9. IX Simulation length and time scales
  10. X Calculation of “strain” in the simulations
  11. XI Fits to simulated time traces
  12. XII Robustness of Simulation results
    1. XII.1 Varying $\sigma_{th}$
    2. XII.2 2$d$ Simulations
    3. XII.3 Next-nearest neighbors
    4. XII.4 Penetration depth effects
  13. XIII Simulated Cluster Properties
  14. XIV Temperature Effects
    1. XIV.1 Transient pump-induced heating
    2. XIV.2 Temperature-dependent time traces
  15. XV Landau theory for SHG near the metal-insulator transition
    1. XV.1 Symmetry-preserving transition
    2. XV.2 Symmetry-allowed components of the SHG susceptibility tensor
    3. XV.3 First approach: direct expansion
    4. XV.4 Second approach: Time-averaged Landau free energy
      1. XV.4.1 SHG in a toy model
      2. XV.4.2 SHG near the metal-insulator transition in CRO

##  I Rotational-Anisotropy SHG fitting

The observed rotational anisotropy (RA) SHG pattern is dominated by the
leading-order electric dipole (ED) contribution to the nonlinear
susceptibility $\chi^{ED}_{ijk}$, defined by:

$P_{i}(2\omega)=\chi^{ED}_{ijk}E_{j}(\omega)E_{k}(\omega),$ (S1)

where $E_{i}(\omega)$ is the incident electric field at frequency $\omega$,
$P_{i}(2\omega)$ is the induced polarization in the sample at frequency
$2\omega$, and a summation over repeated spatial indices is implied. This
contribution to the measured SHG is in general larger than the subleading
magnetic-dipole and electric-quadrupole terms by a factor of the incident
wavelength divided by the lattice spacing Fiebig  _et al._ (2005), which here
is $\approx$500.

Figure S1: RA-SHG data shown in Fig. 1(a) with (a) a fit to Eq. (S4) (this is
the same as the polar plot in Fig. 1(a)) and (b) overlayed with a function
given by (S5) with $\theta_{a}$=$\theta_{c}$= 0.05 rad.

In general, the nonlinear susceptibility tensor describing ED SHG,
$\chi^{ED}_{ijk}\equiv\chi_{ijk}$, is a third-rank tensor with 27 elements.
Requiring the tensor to be invariant under all symmetry operations of the
crystallographic point group $C_{2v}$ ($mm2$), enforcing the permutation
symmetry of indices $j$ and $k$, and taking the axis of inversion symmetry
breaking to be the $b$ axis (see Sec. XV.2) leaves 5 independent non-vanishing
tensor elements:

$\displaystyle\chi_{cbc}=\chi_{ccb},$ (S2)
$\displaystyle\chi_{aba}=\chi_{aab},$ $\displaystyle\chi_{bcc},$
$\displaystyle\chi_{baa},$ $\displaystyle\chi_{bbb}.$

The intensity of the second harmonic is proportional to the norm squared of
the dot product of $\mathbf{P}(2\omega)$ and the unit vector along the
polarization direction of the analyzer $\mathit{\mathbf{e}}(2\omega)$:

$I^{2\omega}\propto|\mathit{\mathbf{e}}(2\omega)\cdot\mathit{\mathbf{P}}(2\omega)|^{2}.$
(S3)

We then calculate the SHG intensity as a function of the polarization angle
$\phi$ of the incident electric field for the cases in which the analyzer is
either parallel or perpendicular to the incident light polarization. In the
normal-incidence geometry employed here, the incident electric field is
restricted to the $ab$ plane, so that all tensor elements with an index $c$ do
not appear in the expression for the SHG intensity. We find:

$\displaystyle
I^{2\omega}_{Par}(\phi)\propto|\chi_{bbb}\mathrm{cos}^{3}\phi+(2\chi_{aba}+\chi_{baa})\mathrm{cos}\phi\mathrm{sin}^{2}\phi|^{2}$
(S4) $\displaystyle
I^{2\omega}_{Perp}(\phi)\propto|\chi_{baa}\mathrm{cos}^{3}\phi+(\chi_{bbb}-2\chi_{aba})\mathrm{cos}\phi\mathrm{sin}^{2}\phi|^{2}.$

These equations give the form of the fit shown in Fig. 1(a), where the tensor
elements are simultaneously fit in the parallel and perpendicular channels.
All time traces were measured with polarization angle $\phi$ = 0 in the
perpendicular channel such that $I^{2\omega}\propto|\chi_{baa}|^{2}$.

###  I.1 Imperfect alignment

The measured RA pattern contains small features that can not be captured by
Equations S4. While these features are not relevant to the discussion in the
main text, for completeness we show that they can be qualitatively captured by
accounting for imperfections in the alignment of the crystal plane to the
probe beam.

The effect of alignment imperfections on the RA-SHG pattern can be
approximated by rotating the $\chi^{ED}$ tensor by small angles before
calculating $I^{2\omega}(\phi)$. This significantly complicates the functional
form of $I^{2\omega}$ and introduces terms containing elements $\chi_{ccb}$
and $\chi_{bcc}$ that are absent when the probe beam is exactly perpendicular
to the $ab$ plane of the crystal.

Figure S1(b) shows the RA data overlayed with a functional form given by Eq.
S5 below. These functional forms account for tilts of the crystal about the
$a$ and $c$ axes away from the ideal orientation in which the $ab$ plane is
exactly normal to the probe beam. Importantly, we find that the bumps in the
perpendicular pattern at 80$\deg$ and 260$\deg$ can be qualitatively accounted
for.

$\displaystyle I^{2\omega}_{Par}(\phi)\propto$
$\displaystyle|(\mathrm{cos}\theta_{c}\mathrm{cos}\theta_{a}\mathrm{cos}\phi-\mathrm{sin}\theta_{c}\mathrm{sin}\phi)(\mathrm{cos}^{2}\phi((2\chi_{aba}+\chi_{baa})\mathrm{cos}^{2}\theta_{a}\mathrm{sin}^{2}\theta_{c}+(2\chi_{cbc}+\chi_{bcc})\mathrm{sin}^{2}\theta_{a})+$
(S5)
$\displaystyle(2\chi_{aba}+\chi_{baa}-\chi_{bbb})\mathrm{cos}\theta_{a}\mathrm{sin}(2\theta_{c})\mathrm{cos}\phi\mathrm{sin}\phi+\chi_{bbb}\mathrm{sin}^{2}\theta_{c}+\mathrm{cos}^{2}\theta_{c}(\chi_{bbb}\mathrm{cos}^{2}\theta_{a}\mathrm{cos}^{2}\phi+(2\chi_{aba}+\chi_{baa})\mathrm{sin}^{2}\phi))|^{2}$
$\displaystyle I^{2\omega}_{Perp}(\phi)\propto$
$\displaystyle|\frac{1}{2}(-\mathrm{cos}\theta_{c}(-2\chi_{aba}-\chi_{baa}-\chi_{bbb}+(2\chi_{aba}+\chi_{baa}-\chi_{bbb})\mathrm{cos}(2\theta_{c}))\mathrm{cos}^{3}\theta_{a}\mathrm{cos}^{2}\phi\mathrm{sin}\phi+$
$\displaystyle\mathrm{cos}^{2}\theta_{a}\mathrm{sin}\theta_{c}\mathrm{cos}\phi((2\chi_{aba}-\chi_{baa}-\chi_{bbb}+(2\chi_{aba}+\chi_{baa}-\chi_{bbb})\mathrm{cos}(2\theta_{c}))\mathrm{cos}^{2}\phi-$
$\displaystyle
2(\chi_{baa}-\chi_{bbb}+(2\chi_{aba}+\chi_{baa}-\chi_{bbb})\mathrm{cos}(2\theta_{c}))\mathrm{sin}^{2}\phi)+\mathrm{cos}\theta_{c}\mathrm{cos}\theta_{a}\mathrm{sin}\phi((-2\chi_{cbc}-\chi_{bcc}-2\chi_{baa}+2\chi_{bbb}+$
$\displaystyle
2(2\chi_{aba}+\chi_{baa}-\chi_{bbb})\mathrm{cos}(2\theta_{c})+(2\chi_{cbc}+\chi_{bcc})\mathrm{cos}(2\theta_{a}))\mathrm{cos}^{2}\phi-$
$\displaystyle(-2\chi_{aba}+\chi_{baa}+\chi_{bbb}+(2\chi_{aba}+\chi_{baa}-\chi_{bbb})\mathrm{cos}(2\theta_{c}))\mathrm{sin}^{2}\phi)+\mathrm{cos}\phi\mathrm{\theta_{y}}(-2\chi_{bcc}\mathrm{cos}^{2}\phi\mathrm{sin}^{2}\theta_{a}-$
$\displaystyle(-2\chi_{cbc}+2\chi_{aba}+\chi_{baa}+\chi_{bbb}+(2\chi_{aba}+\chi_{baa}-\chi_{bbb})\mathrm{cos}(2\theta_{c})+2\chi_{cbc}\mathrm{cos}(2\theta_{a}))\mathrm{sin}^{2}\phi))|^{2}.$

###  I.2 Multiple Domains

If there are multiple domains in the probed region of the sample, the total
SHG intensity at polarization angle $\phi$ is given by

$I^{2\omega}_{tot}(\phi)=\sum_{i}f_{i}I^{2\omega}_{i}(\phi),$ (S6)

where the sum is over unique domains and $f_{i}$ gives the fraction of the
probed region of the sample that constitutes the $i^{\mathrm{th}}$ domain.

By referring to Eq. (S4), we see that if there exist two domains related by
$\phi\rightarrow\phi+$90°, then there exists no angle for which the total SHG
intensity $I^{2\omega}_{tot}(\phi)$ is zero. Domains related by a 90$\degree$
rotation of the polar moment therefore lift the nodes in the RA-SHG pattern
and are not present in the data presented here.

On the other hand, the 180$\degree$ polar domains can only change the
magnitude of $I^{2\omega}$ as a function of position in the $ab$ plane. This
is discussed at length in Ref. Lei  _et al._ (2018). Thus, one cannot reliably
compare the magnitude of the pump-induced drop in $I^{2\omega}$ between data-
sets (it varied between $\sim$4-40% in our measurements, but was typically
$\sim$15%). Nonetheless, at a single spot, the jump in $I^{2\omega}$ is
constant, which indicates that the 180$\degree$ domains are static upon photo-
excitation and that the dominant dynamics are associated with the insulator-
metal transition.

##  II SHG Wavelength Dependence

Figure S2(a) shows the dependence of the RA-SHG pattern on the wavelength of
the probe beam between 820nm and 940nm. Fig. S2(b) shows the value of
$\textrm{max}(I^{2\omega}(\phi))$ in the two channels. The only change in the
pattern comes from the dependence of the tensor elements on the light
wavelength. The same fitting function Eq. S5 can be used to fit the RA-SHG
pattern at each wavelength.

Figure S2: (a) Wavelength dependence of the RA-SHG measurements in both
parallel (left) and perpendicular (right) polarization channels. There is a
90$\degree$ offset between the two channels. (b) The maximal value in the
$I^{2\omega}$ curves from a plotted vs. twice the probe energy.

##  III Time Trace Fits

All time traces are fit to Eq. 1 of the main text. For the data point at the
largest fluence (1.2 mJ/cm2), $\tau$ was restricted to be larger than 100 fs
so that the fitting procedure would converge to a solution. Because the best-
fit value of $\tau$ is found to be $\approx$170 fs (i.e. larger than this
lower bound), we do not expect that this constraint significantly altered the
extracted value. However, it is possible that this constraint resulted in an
underestimation of the uncertainties of the extracted $\tau$ and $\alpha$ at
this fluence.

Figure S3: Attempted fits to the 0.21 mJ/cm2 (top) and 0.28 mJ/cm2 (bottom)
time traces shown in Fig. 3 using various fitting functions. The left column
shows the attempted fit of all time points after $t_{0}$ to a simple
exponential decay with variable decay constant. The second column shows
attempted fits to an error function of variable width with a convolution of
the pump-probe cross correlation. The third column shows the fits to Eq. 1
with $\alpha$ fixed to one. The right column shows the fits to Eq. 1 that are
used in the main text.

There are two timescales in Eq. 1 that correspond to physical processes in the
sample: $\tau$ and a shorter timescale represented by the Heaviside function
that is smaller than, or comparable to our experimental time resolution
$\approx$ 280 fs (c.f. Section V). Fig. S3 shows that the inclusion of only
one of these timescales without the other is insufficient to obtain a good fit
to the measured time traces and thus that the presence of both in Eq. 1 is
necessitated by the data. The left panel of Fig. S3 shows that fits lacking a
short timescale to govern the immediate drop in intensity after the pump pulse
cannot capture the data well. In particular, the fits shown are of the data
after $t_{0}$ to a simple exponential decay with no convolution with the
cross-correlation of the pump and probe pulses $g(w_{0},t)$. Similarly, to
demonstrate that the time constant $\tau$ is necessary to obtain a good fit,
we present attempted fits to Eq. S7 in the second column of Fig. S3. This
function has no exponential decay but includes an error function of variable
width $w$ that approximates the Heaviside function in the limit that
$w\rightarrow$0\. Lastly, the third column displays attempted fits to Eq. 1
with $\alpha$ fixed to one. This demonstrates that, to obtain a satisfactory
fit, the exponential term must have a different weight from the term that has
only the Heaviside function.

$u(t)=1+\\\
\Biggl{[}\frac{1}{2}\left(1+\mathrm{Erf}\left(\frac{2\sqrt{\mathrm{ln}(2)}(t-t_{0})}{w}\right)\right)I_{\infty}\Biggr{]}\ast
g(w_{0},t).$ (S7)

##  IV Recovery time traces

Figure S4: SHG time traces taken over 100s of ps to show the recovery between
pulses.

Figure S4 shows the recovery of $I^{2\omega}$ following the pump pulse. These
traces show that the recovery timescale (hundreds of ps) is much longer than
the timescale for the initial insulator-metal transition, which justifies our
separation of timescales in the percolation model.

##  V Estimates of Time Resolution

In addition to the intrinsic widths of the pump and probe pulses, the time
resolution of the experiment is affected by the angle mismatch between the two
beams. In the limit in which both the pump and probe beams are normally
incident on the sample, the cross-correlation between the two Gaussian
profiles each with a $\sigma$=77 fs (181 fs FWHM) is another Gaussian with
$\sigma$=$\sqrt{2}(77$ fs) = 109 fs (256 fs FWHM). For the geometry employed
in the experiment, we can obtain an upper bound on the time resolution by
adding the time interval between the arrival of the pump pulse at the near and
far ends of the transverse probe profile to the “intrinsic” width of the pump
pulse. This correction gives: $\Delta t$ = FWHMprobesin(15$\degree$)/c = 35
fs. Thus, the obtained upper bound for the FWHM of the cross-correlation
between the pump and probe pulses in the experiment is
$\sqrt{(181)^{2}+(181+35)^{2}}$ fs = 282 fs.

We can also estimate the experimental time resolution from the fitted time
traces and compare it to this upper bound. To do this, we fit the derivative
of the time trace corresponding to the largest fluence (1.2 mJ/cm2) to a
gaussian. If we assume that the drop in SHG intensity occurs much faster than
the time resolution of our experiment, the width of the fitted gaussian is a
measure of our experimental time resolution. (Otherwise, it is an upper
bound). The derivative of the fitted curve for the time trace corresponding to
the largest fluence and the gaussian fit to this derivative are shown in Fig.
S5. The experimental time resolution determined from this procedure is 278 fs
(given by the FWHM of the gaussian). This value is in remarkable agreement
with the upper bound estimated above, 282 fs.

Figure S5: The time-derivative of the fitted curve shown in Fig 3 in the main
text for 1.2 mJ/cm2 (blue points) and a gaussian fit (red) used to determine
the experimental time resolution.

##  VI Calculation of $d_{\gamma}$

We use the following equations to calculate the number of photons absorbed by
a single pulse and the number of unit cells impinged by a single pulse

$N_{\gamma}=(1-R)F\frac{1}{h\nu}\pi(\textrm{hwhm}_{pump})^{2},$ (S8)
$N_{u.c.}=\frac{\pi(\textrm{hwhm}_{pump})^{2}\delta(\omega_{pump})}{abc},$
(S9)

where $R$ is the reflectivity of CRO at 1030 nm, $h\nu$ is the energy of a
single 1030 nm photon, and $\delta_{pump}\approx$ 78 nm is the penetration
depth of the pump laser (c.f. Section VII). We find the number of unit cells
per absorbed pump photon to be $N_{u.c.}/N_{\gamma}$ = 11.5 at the largest
$\tau$.

We can then multiply by the average lattice parameter
$\overline{a}\equiv\frac{a+b+c}{3}$ = 10.1 $\mathrm{\AA}$ to obtain the
average distance between absorbed pump photons $d_{\gamma}$:

$d_{\gamma}\approx\frac{N_{u.c.}}{N_{\gamma}}\overline{a}.$ (S10)

##  VII Penetration depth calculations

To determine the penetration depths of the pump laser (1030 nm) and of the
first- (900 nm) and second- (450 nm) harmonic probe beams, we performed
reflectivity measurements of the sample at each of these wavelengths. To
verify this method, we also measure the reflectivity at 800 nm for comparison
with the penetration depth reported in Ref. Lei  _et al._ (2018). We use these
reflectivity measurements in conjunction with previous measurements Lei  _et
al._ (2018) of the real part of the conductivity, $\sigma_{1}$, to determine
the penetration depths.

The reflectivity $R$ is related to the optical conductivity by the following
expression (in cgs units):

$R=\frac{1+\frac{4\pi}{\omega}(\sigma_{1}^{2}+\sigma_{2}^{2})^{1/2}-(\frac{8\pi}{\omega})^{1/2}\left((\sigma_{1}^{2}+\sigma_{2}^{2})^{1/2}+\sigma_{2}\right)^{1/2}}{1+\frac{4\pi}{\omega}(\sigma_{1}^{2}+\sigma_{2}^{2})^{1/2}+(\frac{8\pi}{\omega})^{1/2}\left((\sigma_{1}^{2}+\sigma_{2}^{2})^{1/2}+\sigma_{2}\right)^{1/2}},$
(S11)

where $\sigma_{1(2)}$ is the real(imaginary) part of the optical conductivity.
From this, we extract $\sigma_{2}(\omega)$ at the desired frequencies, which
we use to compute the penetration depth $\delta$:

$\delta=\frac{c}{(2\pi\omega)^{1/2}}\left((\sigma_{1}^{2}+\sigma_{2}^{2})^{1/2}+\sigma_{2}\right)^{-1/2},$
(S12)

where we have taken $\mu_{1}$=1.

The measured values of the reflectivity $R_{meas}$ are included in Table 1. Of
the relevant wavelengths, only the reflectivity at 1030 nm has been previously
measured Lee  _et al._ (2007), and the value reported here (0.18) is in good
agreement with these prior results. However, for all included wavelengths, the
measured reflectivity is inconsistent with the previously reported values for
$\sigma_{1}$. In particular, using the values for $\sigma_{1}$ in the
literature and the measured reflectivity, Eq. S11 yields a complex value for
$\sigma_{2}$, which is defined to be real. This result is unphysical.
Substituting larger values for the reflectivity, however, yields real values
for $\sigma_{2}$ at all considered wavelengths. To determine the penetration
depth, we use the minimum reflectivity ($R_{det}$) that results in a real-
valued (and thus, physical) $\sigma_{2}$. (In all cases, the difference with
the measured reflectivity was less than or equal to 3$\%$). This procedure is
validated by the penetration depth calculated using this method at 800 nm (73
nm), which is in exact agreement with that reported in Ref. Lei  _et al._
(2018). The calculated penetration depth at 450 nm is also similar to that
reported at 400 nm Lei  _et al._ (2018). The results of the penetration depth
calculations are summarized in Table 1.

$\lambda$ (nm) |  $\sigma_{1}$ (THz) | Rdet(Rmeas) |  $\sigma_{2}$ (THz) |  $\delta$ (nm)
---|---|---|---|---
800 | 681 | 0.19(0.16) | 373 | 73
450 | 1330 | 0.20 | 930 | 31
900 | 614 | 0.19(0.16) | 365 | 76
1030 | 576 | 0.20(0.18) | 576 | 78
Table 1: The penetration depth $\delta$ in CRO at the pump wavelength 1030 nm,
at the first- (900 nm) and second- (450 nm) harmonic probe wavelengths, and at
800 nm for comparison to Ref. Lei  _et al._ (2018). The real and imaginary
parts of the optical conductivity, and the measured and determined values for
the reflectivity that are used to calculate the penetration depth are also
included.

##  VIII Considerations of homogeneous mechanisms

In the main text we have argued for the existence of percolating inhomogeneity
in the PIPT by correlating our measurements with simulations of inhomogeneous
dynamics. While the possibility of a homogeneous PIPT mechanism cannot be
excluded by these results, we argue here that such a scenario is unlikely. The
most essential feature of our data that appears to require the invocation of
spatial inhomogeneity is the non-monotonic fluence dependence of the time
constant $\tau$. If the PIPT is thought to be homogeneous, additional
description is needed to explain this behavior.

One recent study Zong  _et al._ (2019) used a time-dependent Ginzburg-Landau
formalism to describe the non-monotonic fluence dependence of the photo-
induced transition time in a charge density wave compound. In this picture,
the slowdown in $\tau$ appears due to fluctuations with correlation lengths
that diverge near a critical fluence. One might suspect that a similar
description applies here. However, the PIPT observed in CRO appears
substantially different from the previous study in several respects.

First, in the previous study there is only one observed timescale governing
the photoinduced transition due to the coherent displacive motion of the
atoms. It is this timescale that exhibits the slowing down behavior. In
comparison, there exist at least two timescales involved in the transition
process measured in this work (Fig. S3). While the presence of two timescales
is an essential feature of our inhomogeneous model (being associated with the
quasi-instantaneous structural change associated with the absorption of light
and with the subsequent percolation), there is no apparent analogue to the
longer timescale in the homogeneous Ginzburg-Landau picture. Additionally, the
change in the decay time constant observed in the previous work is much
smaller ($\sim$300-400 fs) than that seen here, despite both transitions
involving changes to the crystal lattice. This observation suggests that the
coherent motion of atoms cannot explain the exponential part of the dynamics,
which is consistent with mesoscale structural inhomogeneity in CRO. In
summary, unlike a Ginzburg-Landau description, mesoscopic percolation in CRO
offers a natural explanation for (1) why there exist two timescales in CRO and
(2) why the timescale of the decay in CRO should be significantly longer than
that in the previous study Zong  _et al._ (2019).

We have addressed the possibility of an explanation based on a time-dependent
Ginzburg-Landau theory in detail because it is, to our knowledge, the most
straightforward way in which one might attempt to attribute our observations
to homogeneous dynamics. There are surely other possibilities that we have not
considered. However, we emphasize several reasons why the appearance of
mesoscale inhomogeneity may be expected: (1) the well-documented appearance of
inhomogeneity in the vicinity of first-order insulator-metal transitions, (2)
the central role inhomogeneity plays in the observed colossal
magnetoresistance effect, which is observed in CRO, (3) the observation of
inhomogeneity in the Ti-doped compound McLeod  _et al._ (2021), and (4) that
the nearly binary nature of the $c$ axis change observed in equilibrium is
likely to persist away from equilibrium Perez-Salinas  _et al._ (2022).

##  IX Simulation length and time scales

Although the simulations alone are devoid of any length/time scale, we can
extract crude estimates for the physical scales that correspond to elements of
the simulation by comparing the simulation output with the data. By equating
the peak value of $\tau$ from the measurements with that from the simulation,
we get that one time step in the simulation corresponds to several hundred fs.
Using this value, and equating the percolation speed to the “percolation
speed” of the simulations (1 site per time step), we find that the length
scale corresponding to a single site in the simulation is on the order of 1
nm.

##  X Calculation of “strain” in the simulations

The “strain” in the simulations (that dictates a switch to the $L_{c}$ state
when greater than or equal to $\sigma_{th}$) is calculated according to the
following equations:

$\displaystyle\sigma_{3D}(i,j,k)=$ $\displaystyle
s_{i+1,j,k}+s_{i,j+1,k}+s_{i,j,k+1}+\dots$ (S13) $\displaystyle
s_{i-1,j,k}+s_{i,j-1,k}+s_{i,j,k-1}+\dots$ $\displaystyle\beta($
$\displaystyle s_{i+1,j+1,k}+s_{i-1,j+1,k}+s_{i+1,j-1,k}+\dots$ $\displaystyle
s_{i-1,j-1,k}+s_{i+1,j,k+1}+s_{i-1,j,k+1}+\dots$ $\displaystyle
s_{i+1,j,k-1}+s_{i-1,j,k-1}+s_{i,j+1,k+1}+\dots$ $\displaystyle
s_{i,j-1,k+1}+s_{i,j+1,k-1}+s_{i,j-1,k-1}),$ $\displaystyle\sigma_{2D}(i,j)=$
$\displaystyle s_{i,j+1}+s_{i,j-1}+s_{i+1,j}+s_{i-1,j}+\dots$
$\displaystyle\beta(s_{i+1,j+1}$
$\displaystyle+s_{i-1,j+1}+s_{i-1,j-1}+s_{i-1,j-1}).$

where $\sigma_{2D}$($\sigma_{3D}$) is used for a two-(three-)dimensional
simulated sample. $s_{i,j,k}\in\\{0,1\\}$ is the state of site $(i,j,k)$ where
0 and 1 correspond to the $L_{c}$ and $S_{c}$ states, respectively. These
equations also define the next-nearest neighbor coupling with a coupling
strength given by the parameter $\beta\in$[0,1] (which is set to zero in the
simulations included in the main text).

##  XI Fits to simulated time traces

The fits to the simulated time traces were performed with a procedure almost
identical to the fitting of the experimental data. The only difference is the
removal of the cross correlation factor, which accounts for the temporal
resolution associated with the laser pulses. Fig. S6 shows fits of the
simulated time traces to Eq. 1 of the main text excluding the cross-
correlation term.

Figure S6: Fits (dashed red curves) of simulated time traces shown in the main
text to Eq. 1 without a convolution with the cross correlation factor.

##  XII Robustness of Simulation results

Bootstrap percolation models have been widely studied in the mathematics
literature. In the thermodynamic limit, the entire system becomes activated at
a fraction of initialized sites $p_{c}$ which approaches 0 or 1 depending on
the dimension and the number of nearest neighbors required for activation
Schonmann (1992). However, finite-size effects are found to dominate the
dynamics of even the largest computer simulations. Indeed the threshold of
percolation has been found to approach the thermodynamic limit as slowly as
$\frac{C}{\mathrm{log}\left(\mathrm{log}(n)\right)}$ for a system of size
$n^{3}$ in some of the models that we have studied here Balogh  _et al._
(2012, 2009). Thus, the percolation process undergone during the PIPT studied
here is almost certainly far from the thermodynamic limit, and is best studied
via computer simulation.

In this section we show that variations of the simulation discussed in the
main text produce time traces that preserve the salient qualitative features
of the traces. In particular, the peak in the time constant $\tau$ as a
function of fluence fraction $f$ and its alignment with the saturation of
$I_{\infty}$, as well as the $f$-dependence of $\alpha I_{\infty}$ that is
non-monotonic and goes to zero at extremal values of $f$, are reproduced in
the various simulations shown below. These features should therefore be
understood as rather general consequences of our proposed model as opposed to
artifacts of the specific settings with which the simulation in the main text
was performed.

###  XII.1 Varying $\sigma_{th}$

Figure S7 shows the output of simulations using a 30x30x30 array with
$\sigma_{th}$ changed from 3 as in the main text to 2 and to 4. (No
significant differences were observed between simulations using a 30x30x30 and
a 40x40x40 array.) As may be expected, the value of $f$ at which $\tau$ peaks
and $I_{\infty}$ saturates increases with $\sigma_{th}$.

Figure S7: Simulated time traces and extracted fit parameters with
$\sigma_{th}$ = 4 (top) and $\sigma_{th}$ = 2 (bottom) (with $\beta$ = 0).
Note that the $x$-axes are not the same in the two cases.

While the $\sigma_{th}$ = 4 case retains time traces that appear similar to
those in the data and in the $\sigma_{th}$ = 3 case, the $\sigma_{th}$ = 2
case appears qualitatively different. For the $\sigma_{th}$= 2 case and more
generally for $\sigma_{th}<d$ where $d$ is the dimension of the model, one
expects a sharp transition for large enough systems based on rigorous results
in the mathematics literature Holroyd (2003); Balogh  _et al._ (2012). Our
simulations for higher values of $\sigma_{th}$ seem to indicate the lack of a
sharp transition at least at the system sizes that we have simulated. We take
this as evidence that the experiment is best modeled with $\sigma_{th}$ = $d$.
Nonetheless, in all cases the general behavior of the three fit parameters is
preserved.

###  XII.2 2$d$ Simulations

The broad qualitative features of the simulated time traces are unchanged when
the simulation is limited to two spatial dimensions. The results of the 2d
simulation for $\sigma_{th}$=2 and $\beta$=0 are shown in Fig. S8.

Figure S8: Simulated SHG time traces (a) and extracted fit parameters (b) for
a 2d 30x30 simulation with $\sigma_{th}$=2 and $\beta$=0.

Fig. S8 establishes that our model of percolating inhomogeneity reproduces the
qualitative features observed in the experiment even if the simulation is
confined to two dimensions.

###  XII.3 Next-nearest neighbors

We introduce coupling between next-nearest neighbors by setting $\beta>$ 0 in
Eq. S13. Figure S9 shows the simulated time traces and extracted parameters
for the 30x30x30 simulation with $\sigma_{th}$ = 3 and $\beta$ = 0.1. Again we
find that, while the simulated time traces appear more complicated than the
data, the $f$-dependence of the parameters is qualitatively unchanged.

Figure S9: Simulated time traces and extracted fit parameters for a 30x30x30
array with $\sigma_{th}$ = 3 and $\beta$ = 0.1.

One may be tempted to interpret this as indicative that the next-nearest
neighbor coupling between sites is insignificant in the PIPT in CRO. However,
by increasing both $\sigma_{th}$ and $\beta$, one finds that there are many
combinations of the two parameters (in both $2d$ and $3d$) that appear to
reproduce the qualitative features of the data. The results of several of
these simulations are shown in Figures S10 and S11 for the $3d$ and $2d$ cases
respectively.

Figure S10: Simulated time traces and extracted fit parameters for $3d$
30x30x30 array simulations with (a)(b) $(\sigma_{th},\beta)$=(4,0.1) and
(c)(d) $(\sigma_{th},\beta)$=(6,0.3)  Figure S11: Simulated time traces and
extracted fit parameters for $2d$ 30x30 array simulations with (a)(b)
$(\sigma_{th},\beta)$=(3,0.4) and (c)(d) $(\sigma_{th},\beta)$=(4,0.9)

###  XII.4 Penetration depth effects

In principle, the finite penetration depths of the pump and probe lasers
introduce another type of inhomogeneity into the PIPT. Namely, the number of
(pump and probe) photons absorbed by the sample falls off exponentially as a
function of the depth from the surface. In this section we demonstrate that
the qualitative results of the simulations are preserved when this additional
inhomogeneity is accounted for.

The intensity of the pump beam as a function of depth $z$ falls off
exponentially with the penetration depth $\delta(\omega_{pump})$:

$I_{pump}(z)\propto e^{-z/\delta(\omega_{pump})}.$ (S14)

In the simulations, this effect is captured by scaling the probability of
activating a site at $t=0$ (given by the fluence fraction $f$) by this
exponential factor, where $z$ is now the discrete $z$-coordinate in the
simulated array, and the penetration depth is an unknown number of array sites
(discussed further below).

To account for the finite penetration depth of our SHG probe, we must consider
both the penetration depth of the incident laser at frequency $\omega$ and
that of the generated second harmonic 2$\omega$. The second-harmonic field
generated at a depth $z$ is proportional to the square of the field at
position $z$ with frequency $\omega$, which decays exponentially with the
penetration depth $\delta(\omega)$:

$E^{2\omega}(z)\propto E^{\omega}(z)^{2}\propto e^{-z/\delta(\omega)}.$ (S15)

The generated second harmonic intensity is then partially absorbed by the
sample before returning to the surface. Thus, the light that contributes to
the measured SHG intensity is additionally damped along the $z$ direction due
to the penetration depth of the second harmonic. This leads to an effective
penetration depth of the SHG probe:

$I^{2\omega}_{contr}(z)\propto
e^{-z/\delta(2\omega)}e^{-z/\delta(\omega)}=e^{-z/\delta^{eff}},$ (S16)

where we have defined the effective penetration depth:

$\delta^{eff}\equiv\frac{\delta(\omega)\delta(2\omega)}{\delta(\omega)+\delta(2\omega)}.$
(S17)

To incorporate the effect of this finite-size penetration depth into the
simulations, the simulated SHG intensity, which was previously simply
proportional to the number of $S_{c}$ sites, is now proportional to the number
of $S_{c}$ sites weighted by $e^{-z/\delta^{eff}}$, where $z$ is the discrete
z-coordinate in the array and $\delta^{eff}$ is an unknown number of array
sites (discussed further below).

Using the values for the penetration depth at 900 nm and 450 nm in Table 1, we
obtain an effective probe penetration depth of $\approx$ 22 nm, which is
$\approx$ 3.5 times smaller than the penetration depth of the pump beam,
$\delta(\omega_{pump})$. This ratio being significantly greater than one helps
to explain why simulations neglecting penetration depth effects agree well
with experiment; the probe samples from a much smaller volume of the material
than that which is excited by the pump pulse. Thus, in the region sampled by
the probe, the inhomogeneity due to the penetration depth of the pump is a
small effect.

In incorporating finite penetration depth effects into the simulations, the
ratio $\delta(\omega_{pump})/\delta^{eff}$ must be preserved. However, we do
not know the magnitude of the penetration depths that should be used in the
simulations (in units of number of array sites) to obtain the closest
correspondence to the physical experiment. Therefore, in Figure S12 we present
the results of several simulations in which the magnitude of the penetration
depths is varied while the ratio $\delta(\omega_{pump})/\delta^{eff}$ = 3.5 is
fixed. Simulation output corresponding to several different values of
$\delta(\omega_{pump})$ are shown. We have estimated that each site in the
simulations corresponds to a length on the order of $\sim$1 nm in the
experiment (c.f. Section IX). Thus, in Fig. S12(a) we present the simulation
output with $\delta(\omega_{pump})$=78 sites (and $\delta^{eff}$=22 sites).
Because this value is larger than the size of our simulated array, we further
include the results of simulations with smaller values for the penetration
depth. In part due to the ratio $\delta(\omega_{pump})/\delta^{eff}\gg$ 1, the
effects of the penetration depth on the fluence dependence of $\tau$,
$I_{\infty}$, and $\alpha$ are small.

Figure S12: Simulation output accounting for the penetration depth of the pump
$\delta(\omega_{pump})$ and of the probe $\delta^{eff}$. The ratio
$\delta(\omega_{pump})/\delta^{eff}$=3.5 is fixed in all cases. (a-c) show the
results for $\delta(\omega_{pump})$ set to 78, 9.8, and 2.4 sites,
respectively.

##  XIII Simulated Cluster Properties

Better understanding of the PIPT in CRO and related compounds may come from
future experimental investigations of the photoinduced clusters. In Fig. S13
we present the number of clusters and the average cluster size of the $L_{c}$
and $S_{c}$ states found in our simulations once site-switching has stopped
(using the same simulation settings as used in the main text, i.e. a 40x40x40
array with $\sigma_{th}$ = 3 and $\beta$ = 0). A cluster here is defined as a
collection of sites possessing the same state (i.e. $L_{c}$ or $S_{c}$)
connected by an unbroken chain of nearest neighbor links.

Figure S13: Cluster behavior in the 40x40x40 simulation with $\sigma_{th}$=3
and $\beta$=0. (top) Reproduction of Fig. 3(b) for comparison: $\tau$,
$I_{\infty}$, and $\alpha I_{\infty}$ vs. fluence fraction $f$ (middle)
Average area of the simulated $S_{c}$ and $L_{c}$ clusters at
$t\xrightarrow{}\infty$ vs. fluence fraction (bottom) number of $S_{c}$ and
$L_{c}$ clusters at $t\xrightarrow{}\infty$ in the simulations

Figures S14 and S15 show the time-dependence of the cluster sizes and the
number of clusters, respectively, of both the $L_{c}$ and $S_{c}$ phases in
the simulation presented in the main text (i.e. a 40x40x40 array with
$\sigma_{th}$=3 and $\beta$=0).

Figure S14: The cubed root of the average $L_{c}$ and $S_{c}$ cluster sizes
vs. time steps in the $3d$ simulations with $\sigma_{th}$=3 and $\beta$=0.
Figure S15: The logarithm of the average number of $L_{c}$ and $S_{c}$
clusters vs. time steps in the $3d$ simulations with $\sigma_{th}$=3 and
$\beta$=0.

The sigmoid-like dependence of the cluster size on the fluence, and the time
dependence of the cluster sizes could be observed in future experiments using
time-resolved diffuse x-ray or electron diffraction to verify our model of the
photo-induced insulator-metal transition as a percolation phenomenon.

##  XIV Temperature Effects

To ensure that the sample fully recovered between pulses, we varied the
repetition rate of the laser (both pump and probe pulses) and fixed the probe
pulse to arrive shortly before the pump pulse, while keeping the pump and
probe fluences fixed. If the sample fully recovers between pulses,
$I^{2\omega}$ will exhibit a linear trend with the repetition rate (because
the intensity is proportional to the number of probe pulses), while a
sublinear trend indicates an incomplete recovery of the sample between pump
pulses. Figure S16 shows $I^{2\omega}$ vs. the repetition rate for the largest
pump fluence used in this experiment, 1.3 mJ/$\mathrm{cm}^{2}$. We see that at
repetition rates above the 5 kHz used in the experiment, the sample no longer
recovers between pulses of this fluence. This is verified by performing a
linear fit (with zero offset) to the data points $\leq$ 5 kHz.

Figure S16: $I^{2\omega}$ shortly before the arrival of a 1.2
mJ/$\mathrm{cm}^{2}$ pump pulse vs. repetition rate of the pump and probe
pulses. The red line is a linear fit (with zero offset) to the filled data
points. The data included in the text was taken at 5 kHz.

###  XIV.1 Transient pump-induced heating

From the heat capacity, we can determine an upper bound on the transient
increase in the lattice temperature following the pump pulse. To do this, we
use the following expression:

$q=\int_{T_{i}}^{T_{f}}\mathit{dT}c(T),$ (S18)

where $q$ is the average heat absorbed per unit cell from a single pulse,
$T_{i(f)}$ is the initial(final) lattice temperature, and $c(T)$ is the heat
capacity per unit cell. The heat capacity as a function of temperature is
extracted from Ref. Varadarajan  _et al._ (2007), which used a flux-grown
sample (as was used in our measurements). Specifically, this allows us to
extract coefficients for the $T$ and $T^{3}$ terms in $c(T)$.

The following equation gives the absorbed energy per unit cell $q$, calculated
with a weighted average in which the weight is proportional to the probe
sensitivity:

$\displaystyle q$
$\displaystyle=(1-R)F\frac{V_{uc}}{\delta(\omega_{pump})}\frac{\int_{0}^{\infty}dze^{-z/\delta(\omega_{pump})}e^{-z/\delta^{eff}}}{\int_{0}^{\infty}dze^{-z/\delta(\omega_{pump})}}$
(S19)
$\displaystyle=(1-R)F\frac{V_{uc}}{\delta(\omega_{pump})}\frac{1}{1+\frac{\delta^{eff}}{\delta(\omega_{pump})}},$

where $\delta(\omega)$ = 78 nm is the penetration depth of the pump beam,
$\delta^{eff}$ = 22 nm is the effective penetration depth of the SHG probe
(c.f. Eq. S17), $F$ is the fluence of the pump pulse, and $R$ = .18 is the
reflectivity at the pump wavelength.

The upper bound on the transient temperature $T_{f}$ following the pump pulse
is determined by equating these two expressions for the absorbed energy per
unit cell $q$.

The only remaining parameter in this calculation is the crystal temperature
before the arrival of the pump pulse $T_{i}$. Without laser heating, the
crystal is at a nominal temperature of $\approx$4 K. From Fig. 1(a), we see
that the steady-state temperature increase due to the probe laser is
$\approx$2 K near the transition temperature $T_{MI}$. However, at lower
temperatures the heat capacity is smaller and the steady-state temperature
increase due to the probe beam will be larger. But, because we are interested
in obtaining an upper bound for the transient temperature increase, it
suffices to assume a large initial temperature due to the heating from the
probe laser. Table 2 presents the obtained upper bound on the transient
lattice temperature $T_{f}$ immediately following the pump pulse for various
initial temperatures and for fluences $F_{sat}$ = 0.4 mJ/cm2 and $F_{max}$ =
1.2 mJ/cm2.

$T_{i}$ (K) |  $T_{f}$ @ $F=F_{sat}$ (K) |  $T_{f}$ @ $F=F_{max}$ (K)
---|---|---
10 | 31 | 41
20 | 32 | 42
25 | 34 | 43
Table 2: Upper bounds on the lattice temperature $T_{f}$ immediately following
the arrival of the pump laser pulse for pump fluences $F_{sat}$ = 0.4 mJ/cm2
and $F_{max}$ = 1.2 mJ/cm2 and several initial steady-state temperatures
$T_{i}$.

###  XIV.2 Temperature-dependent time traces

The complete data set corresponding to Fig. 2(b) with pump fluence 0.28
mJ/$\mathrm{cm}^{2}$ is shown in Figure S17(b). Comparison with the data shown
in S17(a) serves to highlight the variation in these measurements. Fig. S17(a)
shows a very long relaxation timescale at temperatures just below the
transition. However, no similar trend appears in Fig. S17(b), where time
traces at temperatures near the transition are flat after several ps and have
values between those at temperatures far above and below the transition. As
mentioned previously, the domain structure along the $c$ axis has an effect on
the measured SHG Lei  _et al._ (2018), and this may be related to the
differences between the datasets presented here.

Figure S17: SHG time traces at varying temperature and the corresponding fit
parameters extracted from fitting to Eq. 1. (a) The pump fluence is fixed at
0.28 mJ/$\mathrm{cm}^{2}$ and the polarization is such that
$I^{2\omega}\propto|\chi^{ED}_{baa}|^{2}$. (b) The pump fluence is 1.2
mJ/$\mathrm{cm}^{2}$ and the polarization is such that
$I^{2\omega}\propto|\chi^{ED}_{bbb}|^{2}$. The two datasets were taken at
different locations on the sample.

##  XV Landau theory for SHG near the metal-insulator transition

We develop a Landau theory for the metal-insulator transition at $T_{MI}$=48
K. (Laser heating is ignored in this discussion.) Our main purpose is to
calculate the in-plane components of the SHG susceptibility tensor,
$\chi^{ED}_{ijk}\equiv\chi_{ijk}$ (where $i,j,k$ take values in $\\{a,b\\}$),
for $T$ near $T_{MI}$.

To describe the transition itself, independent of its relation to SHG, we
postulate that a symmetry-preserving order parameter $\eta$ undergoes a first-
order phase transition at $T=T_{MI}$. In this part of the calculation, we
follow Ref. Kuznetsov  _et al._ (2003), in which a Landau theory is developed
for a symmetry-preserving order parameter (see also Ref. Azzolina  _et al._
(2020)).

The symmetry-preserving order parameter in our calculation is intended to
describe the change in crystal structure at the transition and could be, for
instance, the volume change of the unit cell. The microscopic nature of this
order parameter does not affect the calculations below. Indeed, the
calculation below seems to be equally applicable if the order parameter is
instead taken to describe the change in electronic structure, since this is
also symmetry-preserving. However, as discussed in the main text, experimental
evidence from the Fe-doped compound disfavors either the electronic or
magnetic structural change as explanations for the jump in SHG.

To describe SHG, we provide two phenomenological approaches, neither of which
we attempt to justify microscopically, but both of which yield the same
phenomenological parameterization of the measured quantity $I^{2\omega}$. In
the first approach, we expand the SHG susceptibility in the order parameter
$\eta$. Here we have taken inspiration from the Landau calculation of
elastoresistivity in Ref. Shapiro  _et al._ (2015). In the second approach, we
use a phenomenological time-averaged Landau free energy that includes
frequency dependence, following and extending the approach of Refs. Pershan
(1963); Sa  _et al._ (2000). Ref. Sa  _et al._ (2000) has developed this
approach for the case of a symmetry-breaking order parameter; in this case,
certain SHG tensor components must vanish above the transition temperature by
symmetry, but can be non-vanishing below the transition temperature because
the order parameter reduces the symmetry. In the case we consider, the same
SHG tensor components are allowed both above and below $T_{MI}$ (since the
transition does not break any symmetries), and we describe this in an approach
that goes beyond that of Ref. Sa  _et al._ (2000) in that we explicitly
minimize over frequency-dependent internal degrees of freedom in the time-
averaged Landau free energy.

We organize the calculation as follows. First, we review the symmetry-
preserving transition as developed by Ref. Kuznetsov  _et al._ (2003). Second,
we discuss how symmetries constrain the tensor elements of the SHG
susceptibility. Then we present two phenomenological approaches to obtaining
the jump in $I^{2\omega}$.

###  XV.1 Symmetry-preserving transition

We recall here the Landau theory provided by Kuznetsov et al. Kuznetsov  _et
al._ (2003) for a generic symmetry-preserving phase transition. Taking $\eta$
to be a symmetry-preserving order parameter, one writes

$\mathcal{F}_{0}=a_{1}\eta+a_{2}\eta^{2}+a_{3}\eta^{3}+\eta^{4}.$ (S20)

The cubic term is eliminated by the change of variables $\Delta=\eta+a_{3}/4$,
$\tilde{a}_{1}=a_{1}+a_{3}^{3}/8-a_{2}a_{3}/2$,
$\tilde{a}_{2}=a_{2}-\frac{3}{8}a_{3}^{2}$, yielding
$\mathcal{F}_{0}=\tilde{a}_{1}\Delta+\tilde{a}_{2}\Delta^{2}+\Delta^{4}$,
where a $\Delta$-independent constant has been dropped. The solutions to
$\partial\mathcal{F}_{0}/\partial\Delta=0$ are presented in Ref. Kuznetsov
_et al._ (2003). The final result that we need for the present calculation is
that if $D\equiv 8\tilde{a}_{2}^{3}+27\tilde{a}_{1}^{2}<0$ (note that this
requires $\tilde{a}_{2}<0$), then $F_{3}$ is an asymmetric double well, with
the global minimum being on the left if $\tilde{a_{1}}>0$ and on the right if
$\tilde{a_{1}}<0$. Ref. Kuznetsov  _et al._ (2003) presents three solutions
(in the case $D<0$) called $\Delta_{1}$, $\Delta_{2}$, and $\Delta_{3}$. We
find (contrary to Ref. Kuznetsov  _et al._ (2003)) that the global minimum is
$\Delta_{1}$, i.e., the solution defined by

$\Delta_{1}\equiv-2\rho\cos\frac{\varphi}{3},$ (S21)

where $\rho=\operatorname{sgn}\tilde{a}_{1}\sqrt{|\tilde{a}_{2}|/6}$ and
$\varphi=\arccos\left(\frac{\tilde{a}_{1}}{8\rho^{3}}\right)$ (note that
$\varphi$ is real because of $D<0$).

Figure S18: The global minimum $\Delta_{1}$ for $\tilde{a}_{2}=-1$ (all units
are arbitrary).

To get a first order phase transition for our problem, we fix
$\tilde{a}_{2}<0$ and assume a linear dependence on temperature near $T_{MI}$:

$\tilde{a}_{1}=\tilde{a}_{1,0}(T-T_{MI}),$ (S22)

where $\tilde{a}_{1,0}$ is a temperature-independent constant. Then
$\Delta_{1}$ (and thus also the equilibrium value of the original order
parameter $\eta$) jumps discontinuously at $T=T_{MI}$.

For later convenience, we write the global minimum in an alternate form with
the temperature dependence explicitly included:

$\Delta_{1}(T)=-\sqrt{\frac{|\tilde{a}_{2}|}{2}}\left[\operatorname{sgn}\boldsymbol{(}\tilde{a}_{1,0}(T-T_{MI})\boldsymbol{)}+f\boldsymbol{(}\tilde{a}_{1,0}(T-T_{MI})\boldsymbol{)}\right],$
(S23)

where $f(x)$ is a smooth function defined by

$f(x)=\begin{cases}\frac{2}{\sqrt{3}}\cos\left(\frac{1}{3}\arccos\left[\frac{1}{8}\left(\frac{|\tilde{a}_{2}|}{6}\right)^{-3/2}x\right]\right)-1&0<x<(\frac{2}{3}|\tilde{a}_{2}|)^{3/2},\\\
-f(-x)&-(\frac{2}{3}|\tilde{a}_{2}|)^{3/2}<x<0.\end{cases}$ (S24)

We note in particular the expansion of $\Delta_{1}(T)$ near $T=T_{MI}$:

$\Delta_{1}(T)=-\sqrt{\frac{|\tilde{a}_{2}|}{2}}\left[\operatorname{sgn}\boldsymbol{(}\tilde{a}_{1,0}(T-T_{MI})\boldsymbol{)}+\frac{\tilde{a}_{1,0}(T-T_{MI})}{2\sqrt{2|\tilde{a}_{2}|^{3}}}\right]+O((T-T_{MI})^{2}).$
(S25)

Thus, at the transition temperature, $\Delta_{1}(T)$ jumps discontinuously
with a continuous slope whose sign at $T=T_{MI}$ is the same as the sign of
the discontinuity [i.e., $\Delta_{1}(T)$ resembles either Fig. S18 or the same
figure flipped vertically]. The global minimum $\eta_{0}(T)$ of the original
order parameter is then obtained by the additive shift
$\eta_{0}(T)=\Delta_{1}(T)-a_{3}/4$.

From Eq. (S25), we see that, in order to obtain a consistent expansion to
linear order in $T-T_{MI}$, it may be necessary to include temperature
dependence in some of the coefficients ($a_{2}$ and $a_{3}$) that we took as
constant in Eq. (S20), i.e., set
$\tilde{a}_{2}=\tilde{a}_{2,0}+\tilde{a}_{2,1}(T-T_{MI})+O((T-T_{MI})^{2})$
(with $\tilde{a}_{2}<0$ always) and
$a_{3}=a_{3,0}+a_{3,1}(T-T_{MI})+O((T-T_{MI})^{2})$. Eq. (S25) then would
yield

$\eta_{0}(T)=A+B\operatorname{sgn}(T-T_{MI})+C(T-T_{MI})+D|T-T_{MI}|+O((T-T_{MI})^{2}),$
(S26)

where $A,B,C,D$ are temperature-independent constants [e.g.,
$A=-\sqrt{|\tilde{a}_{2,0}|/2}\operatorname{sgn}\tilde{a}_{1,0}$]. In other
words, $\eta_{0}(T)$ jumps discontinuously and its slope can be different on
either side of the jump.

For the expansion (S20) to make sense, the order parameter must be small
[e.g., we should have $|B|\ll|A|$ in (S26)]. A strongly first order transition
would not necessarily fit into this framework.

###  XV.2 Symmetry-allowed components of the SHG susceptibility tensor

Next, we recall how inversion symmetry along a particular axis is broken in
CRO through the scenario of hybrid improper ferroelectricity Benedek and
Fennie (2011); Harris (2011); Benedek  _et al._ (2012); Nowadnick and Fennie
(2016). It is convenient to consider the components of the polar moment along
the crystallographic axes $\hat{\mathbf{a}}\equiv[110]$ and
$\hat{\mathbf{b}}\equiv[\bar{1}10]$:

$P_{a}\equiv\mathbf{P}\cdot\hat{\mathbf{a}},\
P_{b}\equiv\mathbf{P}\cdot\hat{\mathbf{b}}.$ (S27)

In the following discussion, we use the notation and results of Harris (Ref.
Harris (2011)). While Ref. Harris (2011) focuses on
$\text{Ca}_{3}\text{Mn}_{2}\text{O}_{7}$ and
$\text{Ca}_{3}\text{Ti}_{2}\text{O}_{7}$, we can expect CRO to fit into the
same framework.

The key point is that the full symmetry of the reference tetragonal structure
$I4/mmm$ is broken down to $Cmc2_{1}$ ($\\#36$) by the condensation of an
octahedral tilt mode and an octahedral rotation mode. [We note here that
$Cmc2_{1}$ is an alternative expression for the space group $Bb2_{1}m$
mentioned in the main text; see footnote 13 in Ref. Yoshida  _et al._ (2005).]
Furthermore, there is a symmetry-allowed trilinear coupling between the tilt
mode, rotation mode, and polar moment. There are 8 possible domains: first
there are two choices for the wavevector $\mathbf{q}$ associated with both the
tilt and rotation mode, and then there four choices corresponding to the tilt
mode and rotation mode each have an independently chosen sign. The spontaneous
polar moment lies either along the $a$ axis or $b$ axis with either sign, so
the $8$ domains correspond to four possible directions for the spontaneous
polar moment, i.e., $90$ degree domains. (Recall that we focused in the
experiment on a region of the sample in which any domains present, whether one
or more than one, were aligned along a single axis. However, $90$ degree
domains were observed in other samples.) We note that, up to an overall
constant (which is of no concern since we only wish to establish the axis
along which inversion is broken), we may identify the $x$ and $y$ components
of $\mathbf{P}$ with $Q_{5,1}$ and $Q_{5,2}$ from Ref. Harris (2011). Thus,
again up to an overall constant, our $P_{a}$ and $P_{b}$ are $Z_{5}$ and
$-Y_{5}$ from Ref. Harris (2011). See Eq. (7) of Ref. Harris (2011) for the
trilinear coupling.

The condensation of these order parameters described above occurs at some
temperature much higher than the transition $T_{MI}$ that is our focus. For
our purpose of studying $T$ near $T_{MI}$, we only need to note that inversion
is broken along the axis of the spontaneous polar moment (either $a$ or $b$,
depending on the domain). The symmetries that mix $a$ and $b$ components
together (see Table II and Eq. (2) of Ref. Harris (2011)) are also broken. We
thus see that the symmetry-allowed tensor elements of the SHG susceptibility
are as in (S2).

From here on, we focus on the three tensor elements relevant to the
experimental geometry, namely, $\chi_{aba},\chi_{baa}$, and $\chi_{bbb}$ (note
that $\chi_{aab}=\chi_{aba}$). Also, we consider only the case of a uniform
sample (no domains).

###  XV.3 First approach: direct expansion

For the regime of $T$ near $T_{MI}$ that is our focus, we can obtain a
phenomenological expression for the SHG susceptibility by supposing that it is
analytic in the symmetry-preserving order parameter $\eta$. In particular, we
write

$\chi_{ijk}=A^{(\omega)}_{ijk}+B^{(\omega)}_{ijk}\eta+O(\eta^{2}),$ (S28)

where $(ijk)=(aba),(baa),$ or $(bbb)$ and where
$A^{(\omega)}_{ijk},B^{(\omega)}_{ijk}$ are constants (possibly complex) that
could depend on temperature as well as frequency. This approach has some
similarity to the Landau theory calculation of elastoresitivity in Ref.
Shapiro  _et al._ (2015), though the symmetry considerations here are
different.

It is simplest to assume that $A^{(\omega)}_{ijk},B^{(\omega)}_{ijk}$ are
smooth functions of temperature (at least for $T$ near $T_{MI}$). Then, the
discontinuous jump in the equilibrium value $\eta_{0}(T)$ as $T$ crosses
$T_{MI}$ (see Sec. XV.1) yields a corresponding jump in each $\chi_{ijk}$
(weighted by the appropriate coefficient $B^{(\omega)}_{ijk}$ evaluated at
$T=T_{MI}$). In this way, we have related a first-order, symmetry-preserving
phase transition with a jump in SHG susceptibility.

###  XV.4 Second approach: Time-averaged Landau free energy

We now present another phenomenological calculation of the SHG susceptibilty
for $T$ near $T_{MI}$. Here, the assumption (S28) is replaced by two other
assumptions: (a) neglect of dissipation and (b) the existence of a time-
averaged Landau free energy (see below). This approach yields the same
phenomenological parameterization (of $\chi_{ijl}$) as the approach of Sec.
XV.3 except with real coefficients instead of complex, due to the assumption
(a).

In Ref. Pershan (1963), Pershan proposed a phenomenological “time-averaged
free energy” function that includes frequency dependence and can be used to
calculate nonlinear optical effects. In Ref. Sa  _et al._ (2000), Sa et al.
extended this approach further by including coupling to an order parameter in
the time-averaged free energy. The time-averaged free energy as considered by
Ref. Sa  _et al._ (2000) is a frequency-dependent analogue to the usual
(frequency-independent) equilibrium free energy; it is sufficient for
determining whether or not a particular tensor component of the SHG
susceptibility is or is not constrained to be zero above a symmetry-breaking
transition temperature. Since in our present problem we are considering a
transition that does not break any symmetries, we must develop the approach of
Ref. Sa  _et al._ (2000) further in order to calculate the change in tensor
components that are symmetry-allowed both above and below the transition. To
do this, we use a frequency-dependent analogue to the usual (frequency-
independent) Landau free energy. The essential difference with Sa  _et al._
(2000), then, is that we explicitly minimize over frequency-dependent degrees
of freedom of the system.

####  XV.4.1 SHG in a toy model

As a warm-up, we consider a toy model of SHG (with no temperature dependence
or phase transition): a one-dimensional model without inversion symmetry. The
sample is subject to external electric fields at frequencies $\omega$ and
$2\omega$. We write the amplitudes of the AC electric field and polar moment
at frequency $\omega^{\prime}$ ($\omega^{\prime}=\omega$ or $2\omega$) as
$E^{(w^{\prime})}$ and $P^{(w^{\prime})}$. For a more careful discussion of
the order of limits used in defining these Fourier amplitudes, see Eq. (2.8)
and below in Ref. Pershan (1963).

The AC polar moments, which are the internal variables of the system, appear
in all terms allowed by symmetry. We have already accounted for translation
symmetry by writing the electric fields and polar moments as position-
independent quantities. The only remaining symmetry of this toy model is time
translation. Under time translation by $t^{\prime}$, $P^{(\omega^{\prime})}\to
e^{it^{\prime}\omega^{\prime}}P^{(\omega)}$ and $P^{(\omega^{\prime})*}\to
e^{-it^{\prime}\omega^{\prime}}P^{(\omega)}$, which implies that in any given
term, the sum of frequencies in unconjugated polar moments must equal the sum
of frequencies in conjugated polar moments. For the time-averaged Landau free
energy, we thus write the most general polynomial (up to fourth order) that
satisfies time translation invariance:

$\mathcal{F}_{\text{AC}}^{(\text{toy})}=\sum_{\omega^{\prime}=\omega,2\omega}(a^{(\omega^{\prime})}|P^{(\omega^{\prime})}|^{2}+\frac{1}{2}b^{(\omega^{\prime})}|P^{(\omega^{\prime})}|^{4})+c^{(\omega)}|P^{(\omega)}|^{2}|P^{(2\omega)}|^{2}\\\
+2\text{Re}[d^{(\omega)}P^{(2\omega)*}(P^{(\omega)})^{2}),$ (S29)

in which $a^{(\omega^{\prime})}$, $b^{(\omega^{\prime})}$, $c^{(\omega)}$, and
$d^{(\omega)}$ are constants. Although a complex $d^{(w)}$ is consistent with
symmetry and with $\mathcal{F}_{\text{AC}}^{(\text{toy})}$ being real, an
imaginary part seems to result either in dissipation (contrary to the basic
assumption of the approach Pershan (1963)) or in unphysical effects, so
hereafter we assume that $d^{(\omega)}$ is real. On physical grounds, we also
assume that the parameters in $\mathcal{F}^{(\text{toy})}$ are such that the
AC polar moment vanishes in the absence of an external electric field; for
instance, $a^{(\omega^{\prime})}$, $b^{(\omega^{\prime})}$, and $c^{(\omega)}$
should all be positive.

Including the coupling between the polar moment and the external electric
field, we obtain:

$\mathcal{F}^{(\text{toy})}=\mathcal{F}_{\text{AC}}^{(\text{toy})}+\mathcal{F}_{\text{ext}},$
(S30)

where

$\mathcal{F}_{\text{ext}}=\sum_{\omega^{\prime}=\omega,2\omega}\left(-E^{(\omega^{\prime})}P^{(\omega)^{\prime}*}-E^{(\omega^{\prime})*}P^{(\omega)^{\prime}}\right).$
(S31)

The minimum of $\mathcal{F}^{(\text{toy})}$ (varying the real and complex
parts of $P^{(\omega)}$ and $P^{(2\omega)}$) is then a function of
$(E^{(\omega)},E^{(2\omega)})$ only and coincides with the time-averaged free
energy considered by Ref. Sa  _et al._ (2000).

We proceed to calculate the linear susceptibility and SHG susceptibility in
this toy model. It is convenient to minimize in $P^{(\omega)}$ and
$P^{(2\omega)}$ with $P^{(\omega)*}$ and $P^{(2\omega)*}$ treated as
independent variables. (This is equivalent to minimizing in the real and
complex parts of $P^{(\omega)}$ and $P^{(2\omega)}$.) To present the
calculation more compactly, we define $4$-component arrays:

$\displaystyle\underline{P}$
$\displaystyle=(P^{(\omega)},P^{(\omega)*},P^{(2\omega)},P^{(2\omega)*}),$
(S32a) $\displaystyle\underline{E}$
$\displaystyle=(E^{(\omega)},E^{(\omega)*},E^{(2\omega)},E^{(2\omega)*}).$
(S32b)

The minimization condition may be written as
$\partial\mathcal{F}^{(\text{toy})}/\partial\underline{P}^{*}=0$, i.e.,

$\underline{E}=\frac{\partial\mathcal{F}_{\text{AC}}^{(\text{toy})}}{\partial\underline{P}^{*}}.$
(S33)

The polar moment that minimizes $\mathcal{F}^{(\text{toy})}$ may generally be
expanded

$\underline{P}_{\alpha}=(\chi_{1})_{\alpha\beta}\underline{E}_{\beta}+(\chi_{2})_{\alpha\beta\gamma}\underline{E}_{\beta}\underline{E}_{\gamma}+\dots,$
(S34)

where repeated indices are summed (from $1$ to $4$) and where
$\underline{P}^{s}$ is the spontaneous polar moment and $\chi_{1}$,$\chi_{2}$
are the first- and second-order susceptibility tensors. The SHG susceptibility
is a particular component of $\chi_{2}$:

$\chi_{\text{SHG}}=\frac{1}{2}\frac{\partial^{2}P^{(2\omega)}}{\partial(E^{(\omega)})^{2}}\biggr{\rvert}_{\underline{E}=\underline{0}}=(\chi_{2})_{311}.$
(S35)

To calculate the susceptibilities, we write the $4$-by-$4$ matrix
$\partial\underline{P}_{\alpha}/\partial\underline{E}_{\beta}$ as
$\partial\underline{P}/\partial\underline{E}$ and use a simple matrix
identity:

$\frac{\partial\underline{P}}{\partial\underline{E}}=\left(\frac{\partial\underline{E}}{\partial\underline{P}}\right)^{-1}.$
(S36)

Then we get
$\chi_{1}=(\partial\underline{E}/\partial\underline{P})^{-1}\rvert_{\underline{P}=\underline{0}}=\text{diag}(\chi_{1}^{(\omega)},\chi_{1}^{(\omega)},\chi_{1}^{(2\omega)},\chi_{1}^{(2\omega)})$,
where $\chi_{1}^{(\omega^{\prime})}=1/a^{(\omega^{\prime})}$ [c.f. Eqs. (S29)
and (S33)].

To get the second-order susceptibility, we expand
$\partial\underline{E}/\partial\underline{P}$ to another order:

$\frac{\partial\underline{E}}{\partial\underline{P}}=\chi_{1}^{-1}+\frac{\partial^{2}\underline{E}}{\partial\underline{P}_{\alpha}\partial\underline{P}}\underline{P}_{\alpha}+\dots,$
(S37)

where the omitted terms are quadratic or higher in $\underline{P}$. Then

$\left(\frac{\partial\underline{E}}{\partial\underline{P}}\right)^{-1}=\chi_{1}-\chi_{1}\frac{\partial^{2}\underline{E}}{\partial\underline{P}_{\alpha}\partial\underline{P}}\chi_{1}\underline{P}_{\alpha}+\dots,$
(S38)

which yields the SHG susceptibility:

$\displaystyle\chi_{\text{SHG}}$
$\displaystyle=\frac{1}{2}\left\\{\frac{\partial\underline{P}_{\alpha}}{\partial
E^{(\omega)}}\frac{\partial}{\partial\underline{P}_{\alpha}}\left[\left(\frac{\partial\underline{E}}{\partial\underline{P}}\right)^{-1}\right]_{31}\right\\}\biggr{\rvert}_{\underline{P}=\underline{0}}$
(S39a)
$\displaystyle=\frac{1}{2}\left\\{\chi_{1}^{(w)}\frac{\partial}{\partial
P^{(\omega)}}\left[-\chi_{1}^{(2\omega)}\frac{\partial^{2}E^{(2w)}}{\partial\underline{P}_{\alpha}\partial
P^{(\omega)}}\chi_{1}^{(\omega)}\underline{P}_{\alpha}\right]\right\\}\biggr{\rvert}_{\underline{P}=\underline{0}}$
(S39b)
$\displaystyle=-\frac{1}{2}\chi_{1}^{(2\omega)}(\chi_{1}^{(\omega)})^{2}\frac{\partial^{2}E^{(2w)}}{\partial^{2}P^{(\omega)}}\biggr{\rvert}_{\underline{P}=\underline{0}}$
(S39c)
$\displaystyle=-d^{(\omega)}\chi_{1}^{(2\omega)}(\chi_{1}^{(\omega)})^{2}.$
(S39d)

####  XV.4.2 SHG near the metal-insulator transition in CRO

Next, we combine the above ingredients into a time-averaged Landau free energy
that describes the SHG susceptibility tensor in CRO near the transition at
$T_{MI}$.

As our focus is on the in-plane (i.e., the $a$-$b$ plane) susceptibility, we
ignore the $c$ axis. We write the most general time-averaged Landau free
energy consistent with inversion symmetry along $a$ (but broken inversion
along $b$) and time translation invariance, keeping only up to quartic terms:

$\mathcal{F}=\mathcal{F}_{0}+\mathcal{F}_{\text{AC}}+\mathcal{F}_{\text{coupling}}+\mathcal{F}_{\text{ext}},$
(S40)

where $\mathcal{F}_{0}$ is the Landau free energy for the symmetry-preserving
order parameter $\eta$ [Eq. (S20)], and where

$\displaystyle\mathcal{F}_{\text{AC}}$
$\displaystyle=\sum_{i=a,b}\left[\sum_{\omega^{\prime}=\omega,2\omega}\left(a_{i}^{(\omega^{\prime})}|P_{i}^{(\omega^{\prime})}|^{2}+\frac{1}{2}b_{i}^{(\omega^{\prime})}|P_{i}^{(\omega^{\prime})}|^{4}\right)+c_{ij}^{(\omega)}|P_{i}^{(\omega)}|^{2}|P_{j}^{(2\omega)}|^{2}\right]$
$\displaystyle\qquad+2\text{Re}[2d_{aba}^{(\omega)}P_{a}^{(2\omega)*}P_{a}^{(\omega)}P_{b}^{(\omega)}+d_{baa}^{(\omega)}P_{b}^{(2\omega)*}(P_{a}^{(\omega)})^{2}+d_{bbb}^{(\omega)}P_{b}^{(2\omega)*}(P_{b}^{(\omega)})^{2}],$
(S41a) $\displaystyle\mathcal{F}_{\text{coupling}}$
$\displaystyle=\sum_{i=a,b}\sum_{\omega^{\prime}=\omega,2\omega}\left(\beta_{i}^{(\omega^{\prime})}\eta|P_{i}^{(\omega^{\prime})}|^{2}+\lambda_{i}^{(\omega^{\prime})}\eta^{2}|P_{i}^{(\omega^{\prime})}|^{2}\right)$
$\displaystyle\qquad+2\text{Re}[2\delta_{aba}^{(\omega)}\eta
P_{a}^{(2\omega)*}P_{a}^{(\omega)}P_{b}^{(\omega)}+\delta_{baa}^{(\omega)}\eta
P_{b}^{(2\omega)*}(P_{a}^{(\omega)})^{2}+\delta_{bbb}^{(\omega)}\eta
P_{b}^{(2\omega)*}(P_{b}^{(\omega)})^{2}],$ (S41b)
$\displaystyle\mathcal{F}_{\text{ext}}$
$\displaystyle=\sum_{i=a,b}\sum_{\omega^{\prime}=\omega,2\omega}2\text{Re}[-P_{i}^{(\omega)^{\prime}*}E_{i}^{(\omega^{\prime})}].$
(S41c)

The $a$ and $b$ components of the polar moment and electric field are defined
as in Eq. (S27), now with frequency dependence. The quantities
$a_{i}^{(\omega^{\prime})},b_{i}^{(\omega^{\prime})},c_{ij}^{(\omega)},d_{ijk}^{(\omega)}$
are all temperature-independent constants (as are $a_{2}$ and $a_{3}$ in
$\mathcal{F}_{0}$). The only temperature dependence in $\mathcal{F}$ is in the
term $\tilde{a}_{1}$ that drives the first-order phase transition in $\eta$.
As in the toy model of the previous section, we assume that the constants are
such that the spontaneous AC polar moment is zero.

The linear susceptibility and SHG susceptibilities are found similarly as in
the toy model. The arrays of polar moments and electric fields acquire an
additional index specifying the component in the $a$-$b$ plane:

$\displaystyle\underline{P}$
$\displaystyle=(P_{a}^{(\omega)},P_{b}^{(\omega)},P_{a}^{(\omega)*},P_{b}^{(\omega)*},P_{a}^{(2\omega)},P_{b}^{(2\omega)},P_{a}^{(2\omega)*},P_{b}^{(2\omega)*}),$
(S42a) $\displaystyle\underline{E}$
$\displaystyle=(E_{a}^{(\omega)},E_{b}^{(\omega)},E_{a}^{(\omega)*},E_{b}^{(\omega)*},E_{a}^{(2\omega)},E_{b}^{(2\omega)},E_{a}^{(2\omega)*},E_{b}^{(2\omega)*}).$
(S42b)

The minimization conditions are $\partial\mathcal{F}/\partial\eta=0$ and
$\partial\mathcal{F}/\partial\underline{P}_{\alpha}=0$, the second of which
yields the electric field as a function of polar moment:

$\underline{E}_{\alpha}=\frac{\partial(\mathcal{F}_{\text{AC}}+\mathcal{F}_{\text{coupling}})}{\partial\underline{P}_{\alpha}^{*}}.$
(S43)

We note that the $\eta$-dependent terms in $\mathcal{F}$ may be brought to the
form

$\mathcal{F}_{0}+\mathcal{F}_{\text{coupling}}=a_{1}(\underline{P})\eta+a_{2}(\underline{P})\eta^{2}+a_{3}\eta^{3}+\eta^{4},$
(S44)

where

$\displaystyle a_{1}(\underline{P})$
$\displaystyle=a_{1}+\sum_{i=a,b}\sum_{\omega^{\prime}=\omega,2\omega}\beta_{i}^{(\omega^{\prime})}|P_{i}^{(\omega^{\prime})}|^{2}$
$\displaystyle\qquad+2\text{Re}[\delta_{aba}^{(\omega)}P_{a}^{(2\omega)*}P_{a}^{(\omega)}P_{b}^{(\omega)}+\delta_{baa}^{(\omega)}P_{b}^{(2\omega)*}(P_{a}^{(\omega)})^{2}+\delta_{bbb}^{(\omega)}P_{b}^{(2\omega)*}(P_{b}^{(\omega)})^{2}],$
(S45a) $\displaystyle a_{2}(\underline{P})$
$\displaystyle=a_{2}+\sum_{i=a,b}\sum_{\omega^{\prime}=\omega,2\omega}\lambda_{i}^{(\omega^{\prime})}|P_{i}^{(\omega^{\prime})}|^{2}$
(S45b)

The minimization in $\eta$ thus yields a value $\eta_{0}(T,\underline{P})$
that can be obtained by replacing $a_{1}\to a_{1}(\underline{P})$ and
$a_{2}\to a_{2}(\underline{P})$ in the solution found in Sec. XV.1. Note that
$\eta_{0}(T,\underline{0})=\eta_{0}(T)$ [where $\eta_{0}(T)=\Delta_{1}(T)$
with $\Delta_{1}(T)$ given by Eqs. (S23)-(S25)]. Also note that since
$\underline{P}$ appears at least quadratically in $a_{1}(\underline{P})$ and
$a_{2}(\underline{P})$, we have

$\frac{\partial\eta(\underline{P})}{\partial\underline{P}_{\alpha}}\biggr{\rvert}_{\underline{P}=\underline{0}}=0.$
(S46)

As in the toy model case, we calculate the first-order and SHG
susceptibilities using Eq. (S36) (which is now an $8$-by-$8$ matrix equation).
From
$\chi_{1}=(\partial\underline{P}/\partial\underline{E})\rvert_{\underline{P}=\underline{0}}$,
Eqs. (S41a)-(S41b) and (S43) yield

$\chi_{1}\leftrightarrow\text{diag}(\chi_{1}^{(\omega)},\chi_{1}^{(\omega)},\chi_{1}^{(\omega)},\chi_{1}^{(\omega)},\chi_{1}^{(2\omega)},\chi_{1}^{(2\omega)},\chi_{1}^{(2\omega)},\chi_{1}^{(2\omega)}),$
(S47)

where

$\chi_{1}^{(\omega^{\prime})}=\frac{1}{a^{(\omega^{\prime})}+\beta^{(\omega^{\prime})}\eta_{0}(T)+\lambda^{(\omega^{\prime})}\eta_{0}(T)^{2}}.$
(S48)

By the same manipulations as in the toy model case, we obtain the in-plane
components of the SHG susceptibility tensor:

$\chi_{ijk}=-\left[d_{ijk}^{(\omega)}+\delta_{ijk}^{(\omega)}\eta_{0}(T)\right]\chi_{1}^{(2\omega)}(\chi_{1}^{(\omega)})^{2},$
(S49)

where $(ijk)=(aba),(baa),$ or $(bbb)$.

Dividing out the first-order susceptibilities yields:

$\frac{\chi_{ijk}}{\chi_{1}^{(2\omega)}(\chi_{1}^{(\omega)})^{2}}=A_{ijk}^{(\omega)}+B_{ijk}^{(\omega)}\Delta_{1}(T),$
(S50)

where $A_{ijk}^{(\omega)}$ and $B_{ijk}^{(\omega)}$ are temperature-
independent constants determined by the various constants appearing in
$\mathcal{F}$. The temperature dependence of $\Delta_{1}(T)$ is discussed in
Sec. XV.1.

Furthermore, since in truncating $\mathcal{F}$ to quartic order we have
already assumed $\eta$ to be a small parameter (note that this means
$\tilde{a_{2}}$ must be small), we may expand Eq. (S48) in $\eta_{0}(T)$.
Collecting the various constants, we find that $\chi_{ijk}$ itself (even
without dividing the first-order susceptibilities) also takes the form of the
right-hand side of Eq. (S50) (with different constants $A_{ijk}^{(\omega)}$,
$B_{ijk}^{(\omega)}$).

## References

  * Fiebig  _et al._ (2005) M. Fiebig, V. V. Pavlov, and R. V. Pisarev, J. Opt. Soc. Am. B 1, 96 (2005).
  * Lei  _et al._ (2018) S. Lei, M. Gu, D. Puggioni, G. Stone, J. Peng, J. Ge, Y. Wang, B. Wang, Y. Yuan, K. Wang, Z. Mao, J. M. Rondinelli, and V. Gopalan, Nano Lett. 18, 3088–3095 (2018).
  * Lee  _et al._ (2007) J. S. Lee, S. J. Moon, B. J. Yang, J. Yu, U. Schade, Y. Yoshida, S.-I. Ikeda, and T. W. Noh, Phys. Rev. Lett. 98, 097403 (2007).
  * Zong  _et al._ (2019) A. Zong, P. E. Dolgirev, A. Kogar, E. Ergeçen, M. B. Yilmaz, Y.-Q. Bie, T. Rohwer, I.-C. Tung, J. Straquadine, X. Wang, Y. Yang, X. Shen, R. Li, J. Yang, S. Park, M. C. Hoffmann, B. K. Ofori-Okai, M. E. Kozina, H. Wen, X. Wang, I. R. Fisher, P. Jarillo-Herrero, and N. Gedik, Phys. Rev. Lett. 123, 097601 (2019).
  * McLeod  _et al._ (2021) A. S. McLeod, A. Wieteska, G. Chiriaco, B. Foutty, Y. Wang, Y. Yuan, F. Xue, V. Gopalan, L. Q. Chen, Z. Q. Mao, A. J. Millis, A. N. Pasupathy, and D. N. Basov, npj Quantum Materials 6 (2021).
  * Perez-Salinas  _et al._ (2022) D. Perez-Salinas, A. S. Johnson, D. Prabhakaran, and S. Wall, Nature Communications 13 (2022).
  * Schonmann (1992) R. H. Schonmann, The Annals of Probability 20, 174 (1992).
  * Balogh  _et al._ (2012) J. Balogh, B. Bollobás, H. Duminil-Copin, and R. Morris, Trans Am Math Soc 364, 2667 (2012).
  * Balogh  _et al._ (2009) J. Balogh, B. Bollobás, and R. Morris, Ann. Probab. 37, 1329 (2009).
  * Holroyd (2003) A. E. Holroyd, Probability Theory and Related Fields 125, 195 (2003).
  * Varadarajan  _et al._ (2007) V. Varadarajan, S. Chikara, V. Durairaj, X. Lin, G. Cao, and J. Brill, Solid State Communications 141, 402 (2007).
  * Kuznetsov  _et al._ (2003) A. Y. Kuznetsov, V. P. Dmitriev, O. I. Bandilet, and H.-P. Weber, Physical Review B 68, 064109 (2003).
  * Azzolina  _et al._ (2020) G. Azzolina, R. Bertoni, C. Ecolivet, H. Tokoro, S.-i. Ohkoshi, and E. Collet, Physical Review B 102, 134104 (2020).
  * Shapiro  _et al._ (2015) M. C. Shapiro, P. Hlobil, A. T. Hristov, A. V. Maharaj, and I. R. Fisher, Physical Review B 92, 235147 (2015).
  * Pershan (1963) P. S. Pershan, Physical Review 130, 919 (1963).
  * Sa  _et al._ (2000) D. Sa, R. Valentí, and C. Gros, The European Physical Journal B - Condensed Matter and Complex Systems 14, 301 (2000).
  * Benedek and Fennie (2011) N. A. Benedek and C. J. Fennie, Physical Review Letters 106, 107204 (2011).
  * Harris (2011) A. B. Harris, Physical Review B 84, 064116 (2011).
  * Benedek  _et al._ (2012) N. A. Benedek, A. T. Mulder, and C. J. Fennie, Journal of Solid State Chemistry 195, 11 (2012).
  * Nowadnick and Fennie (2016) E. A. Nowadnick and C. J. Fennie, Physical Review B 94, 104105 (2016).
  * Yoshida  _et al._ (2005) Y. Yoshida, S.-I. Ikeda, H. Matsuhata, N. Shirakawa, C. H. Lee, and S. Katano, Phys. Rev. B 72, 054412 (2005).