# Magnetic ion relaxation time distribution within a quantum well

A. Łopion<EMAIL_ADDRESS>Institute of Experimental Physics,
Faculty of Physics, University of Warsaw, ul. Pasteura 5, 02-093 Warsaw,
Poland  A. Bogucki<EMAIL_ADDRESS>Institute of Experimental
Physics, Faculty of Physics, University of Warsaw, ul. Pasteura 5, 02-093
Warsaw, Poland  W. Kraśnicki  Institute of Experimental Physics, Faculty of
Physics, University of Warsaw, ul. Pasteura 5, 02-093 Warsaw, Poland  K. E.
Połczyńska  Institute of Experimental Physics, Faculty of Physics, University
of Warsaw, ul. Pasteura 5, 02-093 Warsaw, Poland  W. Pacuski  Institute of
Experimental Physics, Faculty of Physics, University of Warsaw, ul. Pasteura
5, 02-093 Warsaw, Poland  T. Kazimierczuk  Institute of Experimental Physics,
Faculty of Physics, University of Warsaw, ul. Pasteura 5, 02-093 Warsaw,
Poland  A. Golnik  Institute of Experimental Physics, Faculty of Physics,
University of Warsaw, ul. Pasteura 5, 02-093 Warsaw, Poland  P. Kossacki
Institute of Experimental Physics, Faculty of Physics, University of Warsaw,
ul. Pasteura 5, 02-093 Warsaw, Poland

###### Abstract

Time-resolved optically detected magnetic resonance (ODMR) is a valuable
technique to study the local deformation of the crystal lattice around
magnetic ion as well as the ion spin relaxation time. Here we utilize
selective Mn-doping to additionally enhance the inherent locality of the ODMR
technique. We present the time-resolved ODMR studies of single
(Cd,Mg)Te/(Cd,Mn)Te quantum wells (QWs) with manganese ions located at
different positions along the growth axis – in the center or on the sides of
the quantum well. We observe that spin-lattice relaxation of Mn2+
significantly depends on the ion-carrier wavefunction overlap at low-magnetic
fields. Interestingly, the effect is clearly observed in spite of very low
carrier density, which suggests the potential for control of the Mn2+ ion
relaxation rate by means of the electric field in future experiments.

ODMR, QW, spin-lattice relaxation, carrier density

##  I INTRODUCTION

The properties of the MBE grown nanostructures are determined by the complex
interplay between the characteristics of used materials and growth conditions.
In particular, the quality of the interfaces between the barrier and QW layer
strongly impacts the optical and electrical quality of nanostructures. The
realistic transition between the barrier and the quantum well material differs
from the idealized image of the perfectly-cut interface both regarding lateral
and the growth direction. The main concern related to the latter issue is the
intermixing, i.e., smoothing of the band profile due to unintentional
migration of the chemical elements across the interface. Even in the case of
negligible intermixing, the interface still can exhibit imperfection in the
lateral direction in the form of a variation of the position of the interface.
The variation of the position may be related to the single monolayer steps
Deveaud  _et al._ (1984) but it can also have a larger extent for rough
interfaces Gaj  _et al._ (1994); Regreny  _et al._ (1987); Grieshaber  _et
al._ (1994); Kossacki  _et al._ (1995).

In the last case, irregularities of the interface introduce severe local
deformation. Therefore magnetic ions near the interface experience different
strain values than ions placed in the center of the QW. The distribution of
the deformation of the crystal lattice at the sites of magnetic ions placed in
various positions in nanostructure can affect the distribution of the ions’
spin-lattice relaxation (SLR) time.

The magnetic ion relaxation can also be affected by interactions with
carriers. The presence of the carrier gas increases the SLR rate Scherbakov
_et al._ (2001a); König  _et al._ (2000), due to effective energy transfer
mediation between the ion and the lattice. In addition, the distribution of
carriers may vary across the nanostructure. Consequently, a non-trivial SLR
time distribution across the nanostructure may also be observed.

All of the above shows that determining the detailed distribution of
properties such as local deformation or distribution of carriers is essential
for better design and fabrication of devices combining electronic and spin
effects. Moreover, this implies that one can utilize the control of the
spatial distribution of the carrier gas density to adjust the properties of
spin-lattice dynamics for spintronics applications.

Studying such properties essentially requires a local approach, while
techniques operating on mean values (such as Electron Paramagnetic Resonance –
EPR) are of limited use. Even the optical techniques, which are usually
considered local, do not exhibit spatial resolution high enough to distinguish
the signal originating from different depths of the QW.

Figure 1:  (a) Schematic picture of the QWs samples used in this work.
Manganese ions in the QW are distributed in various sections along the growth
axis – in either of the two sides of the QW near the interfaces with the
barrier layer (samples UW1498 and UW1499), in the center of the well (UW1500),
and in the entire QW (UW1501); (b) The schematic presentation of the spatial
probability density of the first three heavy hole states in the model sample –
UW1501; (c) Reflectance spectra of studied samples measured at zero magnetic
field.

Here we exploit the fact that by tailoring the profile of the doping by the
magnetic ions, we can restrict their presence to selected parts of the
structure. The hyperfine structure of energy levels of each ion is sensitive
to the deformation of the lattice at its site.

Therefore, the deformation present in this part of the QW might be studied
using the absorption of resonant microwave radiation with an applied magnetic
field Lambe and Kikuchi (1960); Wolos  _et al._ (2004); Bogucki  _et al._
(2022). This absorption can be detected either directly or indirectly by
exploiting the fact that the optical properties of the studied material change
when the paramagnetic resonance occurs. Thus, the resonance can be detected as
a change of optical response, which is the essence of the optically detected
magnetic resonance (ODMR) van Gisbergen  _et al._ (1993); Ivanov  _et al._
(2001, 2007); Tolmachev  _et al._ (2020); Shornikova  _et al._ (2020). The
ODMR technique allows studying observed effects locally in the area probed
optically with micrometer resolution and only in the volume where the magnetic
ions are introduced. The ODMR is particularly suitable for low-dimensional
systems. Spatial selectivity, optical spatial resolution, and high sensitivity
make ODMR a perfect technique for performing measurements on carriers and
excitons coupled to the magnetic ions in nanostructures. The ODMR technique is
especially useful for systems with large exchange interaction between magnetic
ions and photocarriers – a shining example of such a system is (Cd,Mn)Te in
the diluted regime. More generally, the diluted magnetic semiconductors (DMS)
exhibit the giant Zeeman effect Gaj  _et al._ (1994), which connects the shift
of the excitonic energy with the magnetization of the system of magnetic ions.

The transitions between Mn2+ energy levels due to absorption of the microwave
(MW) radiation leads to a decrease in magnetization (evidenced by the decrease
of the giant Zeeman shift), which can also be described in terms of increased
effective spin temperature.

In this work, we combine magnetooptical measurements and ODMR technique to
study interfaces between (Cd,Mg)Te barriers and (Cd,Mn)Te QWs. In contrast to
previous studies van Gisbergen  _et al._ (1993); Ivanov  _et al._ (2001,
2007); Tolmachev  _et al._ (2020); Shornikova  _et al._ (2020), the optical
detection in our ODMR experiment is based on reflectivity measurements, which
allows us to monitor the behavior of not only the ground state exciton but
also higher excited states.

The ODMR technique lets us examine the properties of the structures locally,
in the volume where the magnetic ions are incorporated. This local approach,
in conjunction with the unique design of the samples, is beneficial for
studying the distribution of the deformation and SLR time along the growth
axis of the QW. Restricting the Mn2+ ion incorporation to the specific parts
of the QW allows probing the mentioned properties locally – near the
interfaces and in the center of the well. By exploiting different excitonic
states, we investigate other parts of the structure by varying the overlap
with the magnetic ions. We also use time-resolved ODMR to study spin-lattice
relaxation of the manganese ions incorporated in different positions along the
growth axis.

##  II SAMPLES AND EXPERIMENTAL SETUP

The samples containing single quantum wells were grown by molecular beam
epitaxy on the semi-isolating GaAs substrate with 3.5 $\mu$m CdTe buffer
layer. The 20 nm (Cd,Mn)Te QWs are surrounded by the (Cd,Mg)Te barriers with
magnesium content of about 45%. The barrier underneath the QW is 2 $\mu$m
thick (well above critical thickness), while the top barrier is 50 nm.
Manganese content is equal to about 0.5%. This amount of Mn2+ was chosen to
assure sufficient giant Zeeman effect, but negligible direct ion-ion
interactions and a small bandgap offset (less than 8 meV) Gaj and Kossut
(2010). The samples were designed so that manganese ions are placed at
different positions along the growth axis as shown in Fig. 1(a,b). In order to
verify that all the produced samples actually follow the above layout, we
characterized them using reflectance in a magnetic field. The samples are not
intentionally doped. Nevertheless, the top barrier is only 50 nm thick and due
to the surface acceptors on (Cd,Mg)Te the nonintentional hole density in the
QW is expected Maślana  _et al._ (2003).

For the measurements, the samples were mounted on a holder designed
specifically for ODMR experiments Goryca and Bogucki (2021). The microwave
radiation was provided to the samples with a microstrip antenna. The usage of
an antenna has two main advantages over using a microwave cavity setup. First,
the holder gives easy optical access both for PL and reflectance measurements,
preserving the favored directions between the microwave magnetic and electric
fields, as well as an external magnetic field. Second, it allows us to provide
microwave radiation in a wide range of frequencies from 10 GHz up to about 50
GHz without any specific tuning or changing the cavity. All of the above
allows measuring the paramagnetic resonance in a wide range of magnetic
fields.

During measurements, the samples were immersed in pumped liquid helium
(T$\approx$1.7 K) in a magneto-optical cryostat with superconducting coils
providing the magnetic field values up to 3 T. The whole microwave setup was
tested in a wide range of MW frequencies and amplitudes. The best antenna
performance was selected by measuring the ODMR signal on a reference sample
Bogucki  _et al._ (2022).

We used MW radiation and light illumination in the pulsed mode to avoid
unwanted temperature drifts and other slow setup-related disturbances. The
reflectance spectra were obtained using a filtered supercontinuum laser as a
light source. The impinging light was chopped into pulses with an acousto-
optic modulator (AOM) triggered by the signal of the MW generator (Agilent
E8257D). The relative delay between MW and light pulses was varied in order to
get temporal profiles of the ODMR signal. The width of pulses was a few ms for
MW pulses and tens of $\mu$s for light pulses. The latter determines the
resolution of transients. The pulse timing was controlled with a resolution of
about 10 ns. By keeping constant delays between the pulses, we performed time-
integrated ODMR measurements. The pulsed approach allowed us to correct for
non-resonant thermal drift occurring during long experimental runs. At each
data point, two cases were measured: with the light and MW excitation pulses
in-phase (i.e. overlapping in time) and out-of-phase (i.e. with no overlap in
time). The difference in the intensity of the signal (as defined in Sec.
III.2) between these two situations gave us the ODMR signal robust against
small temperature drifts.

##  III EXPERIMENTAL RESULTS

###  III.1 Magnetooptical measurements

Figure 2: (a)Reflectance spectrum measured for different magnetic fields for
sample UW1498; E1HHi denotes excitonic states, $\sigma^{+}/\sigma^{-}$ denotes
detection polarization; (b) Giant Zeeman splitting of the different excitonic
complexes observed for sample UW1498 vs. magnetic field

The magnesium content of the barrier layers was verified using standard
reflectance measurements. We observed barrier exciton at around 2.6 eV, which
corresponds to magnesium content equal to 46%, as
$E_{\mathrm{g}}^{\mathrm{Cd}_{1-x}\mathrm{Mg}_{x}\mathrm{Te}}(x)=Bx+E_{\mathrm{g}}^{\mathrm{CdTe}}=(1.85x+1.606)[\mathrm{eV}]$
Waag  _et al._ (1993). The value of manganese content was determined by
fitting the modified Brillouin function Gaj and Kossut (2010); Gaj  _et al._
(1994, 1979) with the parameters from ref. Gaj  _et al._ , 1994 to the exciton
giant Zeeman shift obtained for the reference sample UW1501 (the Mn2+ ions
present in the whole QW). The Mn2+ content was confirmed to be equal to 0.5%.

The reflectance spectra for each studied QW exhibit several distinct features
related to excitonic complexes of electrons and holes at different QW sub-
bands, see Fig.1(c). The most pronounced feature at the lowest energy, formed
by the lowest energy electron and heavy-hole subbands, corresponds to the
excitonic state denoted as E1HH1. The energy of this exciton is determined by
several factors, such as the width of the QW (the same for all the samples),
but also the overlap between the excitonic wavefunction and incorporated Mn2+
ions, since the bandgap of (Cd,Mn)Te is larger than that of pure CdTe Gaj and
Kossut (2010). For samples with magnetic ions placed only on the sides of the
QW, the overlap between exciton wavefunction and magnetic ions gives small
addition to ground state energy. In contrast, in the case of the sample with
ions placed in the center of the QW, the overlap is significantly larger and
results in a respectively larger increase of the ground state energy. Finally,
manganese ions present in the whole well raise the QW transition energy by
about 8 meV, which is in agreement with 0.5% Mn content determined from the
giant Zeeman effect Gaj and Kossut (2010). Thus, excitonic states in all the
samples are shifted towards higher energies in comparison to the QW without
any Mn2+ ions. For all samples, features corresponding to a number of excited
states (E1HH$i$, $i=2-4$) are less pronounced but still visible. It is worth
noting that we can observe the transition E1HH2, which should be forbidden in
the perfectly symmetric QW. Its visibility suggests asymmetry of the wells’
potential Cibert  _et al._ (1993); Vanelle  _et al._ (1996), observed even for
the samples UW1500 and UW1501, which were designed with a symmetric Mn-doping
profile.

Figure 2(a) shows an example map of the reflectance spectrum measured in the
magnetic field applied in the Faraday configuration. Features corresponding to
excitonic complexes – ground state neutral exciton and a number of excited
states – are well-pronounced. The observed variety of the excitonic
transitions exhibits different behavior in the magnetic field. The samples
with the magnetic ions present on sides near the interfaces have a higher
wavefunction-magnetic-ion overlap for higher states. This leads to a stronger
coupling to the source of magnetization than for the ground state, causing a
more pronounced Zeeman splitting, see Fig. 2(b). This result shows qualitative
agreement with the model and further confirms that the manufactured samples
are structured in accordance with our design. However, for the states higher
than HH2, behavior in the magnetic field becomes more complicated as heavy
holes are mixed with the light holes. The first light hole state is expected
to be about 20 meV above the HH1 state, due to the axial strain induced
splitting Bogucki  _et al._ (2022) and different confinement energy. This
splitting is caused by the deformation $\varepsilon_{\parallel}$ of the QW
created by lattice mismatch with barrier Bogucki  _et al._ (2022):
$\Delta_{LHHH}=2b(1+\frac{2C_{12}}{C_{11}})\varepsilon_{\parallel}+\mathrm{const}$,
where $b=-0.94$ eV is shear deformation potential, $C_{11},C_{12}$ are elastic
stiffness constants and $\frac{2C_{12}}{C_{11}}|_{1.6K}=1.39886$

Additionally, figure 1(b) shows that samples UW1498 and UW1499 with Mn2+
incorporated at two opposite interfaces have slightly different energies of
E1HH1 exciton. They should be symmetric and have the same optical properties
in the ideal case. However, the energy position of excitonic states differs
here, which suggests a different overlap between carriers and Mn2+ volume.
Those differences are relatively small but clearly visible and become even
more pronounced in the magnetic field. The Zeeman splitting of the exciton
ground state in the sample UW1499 (with manganese ions near the top interface)
is almost twice as high as splitting in sample UW1498 (bottom interface) see
Fig.3(a). Zeeman splittings for the UW1500 and UW1501 are even higher due to
the higher wavefunction overlaps.

Figure 3: (a) Splitting in the magnetic field of the E1HH1 state measured for
different samples at 1.7 K (b) An overlap between the excitonic ground state
and the magnetic ions volume obtained from the measurement of Zeeman shift vs.
values calculated with consideration of electrostatic potential caused by the
presence of the hole gas of density equal to $0.13\times 10^{11}\,cm^{-2}$
with the spatial distribution given by ground heavy hole state probability.

Different overlap between magnetic ions volume and excitonic states derived
both from different Zeeman splittings and zero field energies, along with the
presence of E1HH2 transition, suggest the asymmetry of the QWs potential. At
least two effects can explain the discrepancy between the Zeeman splittings in
nominally symmetric samples: asymmetry of the top and bottom interfaces or
built-in electric field.

The first one corresponds to the non-abrupt, intermixed character of the
interfaces. The formation process of interfaces depends on material properties
and growth conditions. The QW asymmetry is caused by one of the intermixing
mechanisms – segregation: the growing material is dragged towards the outer
layers, along the growth axis Grieshaber  _et al._ (1996); Gaj  _et al._
(1994). Thus, magnesium from barrier material is mixed into the QW material on
the first interface and pushed out on the second one. This shift of
composition moves potential in the first part of the well upwards, shifting
the carrier wavefunctions away from the first interface.

The second effect corresponds to unintentional p-type doping. Presented QWs
are located near the surface of the sample (the top barrier is 50 nm thick).
In such a case, the surface states act as acceptors, causing the presence of
hole gas in QW Maślana  _et al._ (2003). However, in the reflectance spectrum,
we do not observe any suggestion of charged excitonic features below the
neutral exciton ground state energy. From that fact we can conclude that the
carrier density in the QW, for studied samples is low – i.e. estimated below
1011 cm-2 Kossacki  _et al._ (1999); Kossacki (2003); Kossacki  _et al._
(2004) – yet still significant. The positive sign of charge carriers was
confirmed using the spin-singlet triplet transition in low magnetic fields as
in Łopion  _et al._ (2020); Kossacki  _et al._ (2004) (not shown).

The presence of carrier gas in the QW (holes) and negatively charged states at
the sample’s surface result in an electrostatic field in the QW and the top
barrier. This field modifies wavefunctions of carriers confined in QW. The
modification is most pronounced for heavy holes due to their high effective
mass. The hole gas builds a potential that drags holes wavefunctions towards
the top barrier and electrons towards the bottom barrier. Consequently,
wavefunction overlap with the manganese volume is higher in the sample with
ions located near the top barrier, where the modified potential is deeper
(UW1499).

We analyzed the plausibility of the above scenarios using numerical
simulations. Therefore, nominal QW potential was calculated from composition
using the same parameters as in Ref. Bogucki  _et al._ , 2022. Moreover,
modification of the potential was added based on the intermixed interfaces or
on the electric field, which is rooted in carrier distribution along the
z-axis of the sample (holes in the QW and electrons on the surface). The
calculation of the wavefunctions perturbated by the electric field were
performed with the use of a single-step Schrodinger-Poisson Solver, which is a
sufficient approach for low carrier densities. We calculated overlaps between
wavefunctions and the magnetic ion volumes for each sample as integrals over
the proper part of the QW (A,B,C or the whole QW) of the calculated
wavefunctions density.

We found that the change of the potential corresponding to the intermixed
interfaces alone is insufficient to explain observed differences in the Zeeman
splitting when taking reasonable intermixing coefficients below 1 nm Gaj  _et
al._ (1994). At the same time, calculations for the scenario with the electric
field involved show excellent agreement with the measured data when assuming
carrier density equal to 0.13$\times$1011 cm-2. In the first approximation we
assumed that carrier density is the same for all the samples. That assumption
seems to be reasonable as the design of the samples is very similar –
especially they have the same thickness and composition of the top barrier
layer but generally the carrier density can slightly vary between the samples.
Nontheless from the very good agreement for the samples UW1498 and UW1499 we
conclude that for both of the samples the carrier density is very close to the
one used in the calculation. The noteworthy is the fact that the Zeeman
splitting obtained for the sample UW1500 (the one with the manganese ions
placed in the middle) is slightly higher than it should be in the ideal,
symmetric case. That cannot be obtained with the change of the carrier density
as the addition of built-in electric field can only decrease the overlap
between the ground state wavefunctions and the Mn ions volume in the centre of
the QW.

The calculated overlap was defined as integral of state probability over the
part of the QW where Mn2+ ions are present. Assuming 100% overlap between HH1
wavefunction and magnetic layer for sample UW1501, which causes a Zeeman
splitting equal to 13.2 meV in 3 T, we can calculate the overlap between HH1
and magnetic layer for the other samples. This overlap is obtained as a ratio
of experimentally determined Zeeman splitting for the analyzed sample and the
sample UW1501, see Fig. 3(b). We can compare those values with the HH1 – Mn2+
overlap. As for sample UW1501, where manganese ions are present in the entire
well, this overlap is equal to almost 1. For other samples, the overlap is
respectively lower. Finally, we conclude that observed differences of giant
Zeeman shift between samples UW1499 and UW1498 correspond to the electric
potential built in the samples.

###  III.2 ODMR measurements

Figure 4: (a) Reflectance spectrum measured for the different magnetic fields
in the presence of microwave radiation. In the resonant magnetic field (here
0.57 T for 15.85 GHz) the energy shift is observed for all excitonic complexes
– (e.g., see dashed rectangle); (b) Reflectance spectrum collected in magnetic
field 0.57 T with resonant microwave radiation (15.85 GHz) with pulses in
positions ON and OFF; insets represent two different excitonic states – ground
state E1HH1 and one of the excited states E1HH3; (c) Normalized ODMR signal
measured for all visible excitonic complexes of UW1498 sample.

Figure 4(a) shows the reflectance spectra measured versus the magnetic field
under microwave radiation. All the excitonic lines in the reflectance spectrum
exhibit sensitivity to the applied microwave radiation for magnetic resonance
fields (near $\pm$0.5 T). We observe that in the magnetic resonance fields,
all the features are shifted towards the energies corresponding to lower
absolute values of the magnetic field, see a dashed rectangle in Fig. 4(a).
Furthermore, the shift of the excited states under microwave radiation is more
pronounced than the shift of the ground state, see Fig. 4(b). That effect
corresponds to a larger change of the Zeeman shift for excited complexes.

We define the ODMR signal as the difference in the energy position of the
excitonic line in the reflectance spectrum for microwave and light pulses
overlapping (ON) and shifted (OFF). An example of a representative ODMR signal
as a function of magnetic field for fixed microwave frequency is presented in
fig. 4(c). The shape of the signal versus magnetic field obtained for all the
excitonic states visible in the single QW is the same, see Fig. 5(a). This
fact clearly shows that all the complexes are probing the ensemble of the
magnetic ions with the same properties.

Figure 5: (a) Normalized ODMR signal vs. magnetic field extracted from E1HH1
state for all measured samples (b) ODMR signal vs. magnetic field measured for
the UW1498 QW with the simulated ODMR signal. The red spiked curve shows
calculated magnetic transitions without any line broadening. The blue curve is
a result of the convolution of the red-lined spectrum with the Gaussian
kernel. (c) Compilation of the simulated shapes of ODMR signal as a function
of magnetic field for different values of the spin Hamiltonian D parameter.

We also observe that the shape of the ODMR signal obtained for different
samples is similar, see Fig. 5(a). The main factor that determines the overall
shape of the ODMR signal as a function of the magnetic field at low
temperatures is the deformation of the crystal lattice surrounding probed
magnetic ions, Bogucki  _et al._ (2022). The lattice mismatch between the
layers of the structures results in nonzero strain-induced axial-symmetry spin
parameter D of the Hamiltonian of the Mn2+ ion Qazzaz  _et al._ (1995).
Angular resolved OMDR measurements can determine the value of the spin
Hamiltonian parameter D. In the previous work Bogucki  _et al._ , 2022 we
assumed homogenous strain within the QW, and we probed whole of its volume
with preferential sensitivity of the central region. In this work we test this
assumption by studying QW with Mn2+ restricted to the specific parts along
axis growth, thus with preferential sensitivity in different parts of the QW.

As the bottom barrier layer is 2 $\mu$m thick we assume that it is relaxed on
the interface with the QW and it governs the lattice mismatch for the QW
layer. The nominal deformation coming from this lattice mismatch is
$\varepsilon_{\parallel}=\frac{a_{barrier}}{a_{QW}}-1$ where $a_{barrier}$ and
$a_{QW}$ are lattice constants for barrier and QW layers, respectively. In the
case of (Cd,Mg)Te/(Cd,Mn)Te QW with Mg content about 46% and Mn content about
0.5% the lattice constant are $a_{barrier}=6.45248$ Å and $a_{QW}=6.48028$ Å
deformation is $\varepsilon^{300K}_{\parallel}=-4.195\,\permil$. Taking into
account the different thermal expansion coefficients of CdTe and substrate
GaAs the deformation for measurement temperature is
$\varepsilon^{1.6K}_{\parallel}=-4.837\,\permil$. As it was presented in
Bogucki  _et al._ , 2022, following the equation
$D=-\frac{3}{2}G_{11}\left(1+\frac{2C_{12}}{C_{11}}\right)\varepsilon_{\parallel}$
such value corresponds to the D parameter about 1250 neV. All the constants
are the same values as presented by Bogucki et al. in Ref. Bogucki  _et al._ ,
2022.

We compare measured ODMR spectra with the numerically obtained positions and
intensities of absorption lines based on the Mn2+ spin Hamiltonian like in
ref. Bogucki  _et al._ (2022). The finite linewidth is obtained using the
phenomenological broadening approach. In figure 5(b) experimental ODMR signal
vs. magnetic field measured for the sample UW1498 is shown along with the
calculated one. Each of the 30 transition lines (marked in red) is convoluted
with a gaussian kernel of FWHM= 34 mT (blue). We satisfactorily reproduce the
measured ODMR spectra for a D parameter value of 1250 neV which is in a good
agreement with the strain nominally present in the structure.

Mention worthy is the sensitivity of the shape of the ODMR signal vs. magnetic
field to changes of the D parameter as is presented in figure 5(c). The
agreement with the D corresponding to the nominal lattice mismatch suggests
that the strain is distributed equally along the growth axis of the QW layer
for each of the samples.

Figure 6: Normalized time-resolved ODMR signal for E1HH1 and E1HH3 excitonic
states for samples: (a) UW1500 and (b) UW1499 for resonant frequency f=10.15
GHz. The observed ODMR amplitude relaxation time is the same regardless of the
excitonic state used for the probe. Figure 7: Time-resolved ODMR signal
amplitude for samples UW1499 and UW1500 extracted from E1HH1 state for
different resonance frequencies: (a) 10.15 GHz and (b) 41.9 GHz.

###  III.3 Time-resolved ODMR measurements

We have performed time-resolved ODMR measurements on each of the presented
samples and the visible excitonic complexes. The spin-lattice relaxation (SLR)
time measured on the higher excitonic states has the same value as those
obtained on the ground state. The example of the temporal profiles is shown in
Fig. 6(a,b). The SLR time measured on the excited states is the same as for
the ground state for all the samples.

As expected, the SLR time depends on the magnetic field (resonant MW
frequency), see Fig. 7(a,b). We observe a significant shortening of the SLR
time for higher magnetic fields, where the differences between samples become
imperceptible.

The analysis of measured SLR times as a function of magnetic field for
different samples reveals non-trivial effects. In figure 7(a), temporal
profiles measured for 10.15 GHz resonance for samples UW1500 and UW1499 are
compared. This frequency corresponds to the paramagnetic resonance in a
magnetic field equal to approximately 0.36 T. In this case, we observe a
significant difference in the SLR times. The situation is noticeably different
when it comes to measurements in higher magnetic fields – temporal profiles
obtained for 41.9 GHz resonance (about 1.5 T) are shown in Fig. 7(b). Here,
the SLR times are the same for both samples.

A summary of the measured SLR times for all samples and the magnetic fields is
presented in figure 8(a). For all analyzed samples, a maximum value of SLR
time is observed around 0.6 T. The SLR times decrease with the magnetic field
and follow the same curve for all the samples for higher magnetic fields.
However, the measurements at lower magnetic fields show notable differences
for different QWs. Especially the samples with Mn2+ ions placed in the center
of the QW (UW1500 and UW1501) have significantly shorter SLR times than the
samples with ions placed on the QW interfaces. This finding suggests an
essential role of the Mn2+ – HH1 probability density overlap in the relaxation
mechanisms at lower magnetic fields.

##  IV Discussion of spin relaxation

Figure 8: (a)SLR time versus the magnetic field. Black arrows mark magnetic
fields for which temporal profiles are presented in Fig.7(a) and (b); shadowed
area in the background is a guide for the eye. As the magnetic field
increases, variation of the SLR times between the samples decreases; (b) mean
HH1 – Mn2+ overlap per manganese volume vs. SLR time for $B=0.36$ T

What we observe is that SLR depends on the applied magnetic field following
the curve shown in figure 8(a). The observed change of behavior of the SLR
time in the magnetic field suggests at least two regimes with different
leading mechanisms causing SLR relaxation.

In general, spin dynamics in DMS are controlled by the energy and spin
transfer between the Mn2+ system, carriers, and the crystal lattice. In our
experiment, we heat up the Mn2+ system resonantly with the MW radiation. The
Mn2+ system can relax directly – with spin and energy dissipated to the phonon
bath or indirectly via carriers present in the structure. Measured SLR is a
result of those two paths. The rates of these energy transfers depends on the
effective temperatures of the systems and densities of states available for
scattering.

$\frac{1}{\tau_{\mathrm{SL}}}=\frac{1}{\tau^{0}_{\mathrm{SL}}}+\frac{1}{\tau_{h-Mn}}$
(1)

We denote here the base SLR time of Mn ions spins and energy directly to the
lattice as $\tau^{0}_{\mathrm{SL}}$. In the regime of the high magnetic
fields, this mechanism takes the leading role as its rate is determined by the
increasing density of phonon states allowed for scattering at the energy
corresponding to the Zeeman splitting. It is well established that the
magnetic field dependence of this mechanism is power of B:
$1/\tau^{0}_{\mathrm{SL}}\sim B^{n}$ Strutz  _et al._ (1992, 1993); Goryca
_et al._ (2015). The exponent of B depends on the model applied for
calculating the density of phonon states and the type of scattering mechanism
(direct or indirect). Such increase of the SLR rate of Mn2+ with the magnetic
field was studied for manganese ions in different systems: bulk crystals,
quantum wells, and quantum dots Scherbakov  _et al._ (2001a); Strutz  _et al._
(1992); Goryca  _et al._ (2015); Strutz  _et al._ (1993). In our particular
case, for $B>1$ T we find that the relaxation time is shorter upon increasing
the magnetic field what is in the agreement with previous reports. In this
regime all studied samples exhibit the same SLR time, which is expected due to
the same phonon spectrum.

Upon approaching the low field limit, the phonon based mechanism is less and
less efficient, and is surpassed by another mechanism: spin and energy
scattering via the carrier’s system. The ultimate intrinsic limitation for the
capacity of this channel is given by the carrier spin lifetime
$\tau_{\mathrm{S}}$, which is very short Scherbakov  _et al._ (2001a); Akimov
_et al._ (2006). As the $\tau_{\mathrm{S}}$ is very short, in practical
situations this channel is limited by the carrier gas-Mn system scattering
time $\tau_{\mathrm{h-Mn}}$.

Conceptually, the simplest mechanism of energy transfer between Mn and carrier
system is the spin flip-flop König  _et al._ (2000); Scherbakov  _et al._
(2001a); Akimov  _et al._ (2006), in which a Mn2+ ion and a carrier undergo a
simultaneous spin flips. The energy conservation requires that the two
involved states of the carrier are separated by the energy equal to the Zeeman
splitting of the Mn2+ ion. Since the initial and final carrier states
correspond to different spin subbands, the availability of the populated
initial states crucially depends on the Fermi level $\varepsilon_{\mathrm{F}}$
(and thus the carrier density). The giant Zeeman splitting of the valence band
is equal $E_{\mathrm{z,c}}$ =
$\mu_{\mathrm{B}}g^{\mathrm{eff}}_{\mathrm{c}}B$, with effective g-factor
$g^{\mathrm{eff}}_{\mathrm{c}}$, which is extraordinarily big – i.e. for
sample UW1501 $g^{\mathrm{eff}}_{\mathrm{h}}\sim 110$. Due to the large value
of the $g^{\mathrm{eff}}_{\mathrm{c}}$, the splitting $E_{\mathrm{z,c}}$
exceeds $\varepsilon_{\mathrm{F}}$ even in relatively small magnetic fields.
In such a case, only the lower spin subband is occupied, while the upper band
is empty. As the Zeeman splitting of Mn2+ ions ($E_{\mathrm{Z,Mn}}$) is
relatively small (g-factor $=2$), the energy transfer from magnetic ions to
the carriers becomes impossible, and the $\tau_{\mathrm{h-Mn}}$ rate decrease
with the increase of the magnetic field. Different Zeeman splittings of Mn2+
and carriers leads to decoupling of both systems at higher magnetic fields.
The maximum of the SLR via carrier channel is in zero magnetic field, while
the separation between energy levels of Mn2+ ions and spin subbands are
similar. For $B>0$, the SLR rate via this channel drops due to the detuning of
the energy levels of those two systems. As the carriers sub-system has a
finite effective temperature also the states with higher k-vectors are
occupied. Those states can take part in the relaxation even if the separation
of sub-bands in the magnetic field exceeds $E_{z,Mn}$ in the magnetic field. A
very similar situation is when the carrier density is sufficiently large, and
the Fermi level exceeds the $E_{\mathrm{Z,c}}$. The flip-flop interactions
model leads to the $\tau_{\mathrm{h-Mn}}$ dependence on the magnetic field
$B$, carrier density via Fermi level $\varepsilon_{\mathrm{F}}$, carriers’
spins scattering time $\tau_{\mathrm{S}}$ and effective temperatures for both
subsystems – magnetic ions $\beta_{\mathrm{Mn}}$ and carriers
$\beta_{\mathrm{c}}$ Scherbakov  _et al._ (2001a); König  _et al._ (2000):
$\tau_{\mathrm{h-Mn}}=\tau(B,\varepsilon_{\mathrm{F}},\tau_{\mathrm{S}},\beta_{\mathrm{c}},\beta_{\mathrm{Mn}})$.

Aforementioned mechanism is present both in p- and n-doped systems. In
comparison, the holes have a higher effective mass $m_{\mathrm{eff}}$ than
electrons, causing Fermi level for the density of holes $n$ to be lower than
for the same density of electrons $n$
($\varepsilon_{\mathrm{F}}=\frac{n\pi\hbar^{2}}{m_{\mathrm{eff}}}$).
Nevertheless, the interactions with holes can provide more efficient
relaxation channel Akimov  _et al._ (2006); Scherbakov  _et al._ (2001b).
Holes have a stronger effect on the SLR as they have a larger density of
states (due to the larger effective mass), and (p-d) exchange interaction of
Mn2+ ions with holes is stronger than with electrons Gaj and Kossut (2010).
Moreover, there may also be a second process that allows SLR of Mn2+ ions via
hole gas. It was considered by Akimov et al. Akimov  _et al._ (2006) that
while for the electrons there are only flip-flop transitions, for the holes,
the transfer of spin may be possible during the relaxation in a single Zeeman
subband due to the mixing of heavy and light hole states, proportional to
$\sim k^{2}$. This mechanism might be an even more effective relaxation
channel than flip-flop transitions if Fermi level is too low for the latter
one.

Summing up: observed drop of the SLR time in low magnetic fields may be
explained by the interactions with hole gas. Moreover, for the magnetic field
below 0.6 T, we observe the differences in SLR times between the studied
samples. Those differences are correlated with the mean overlap between the
manganese ions and the local wavefunction density of HH1, which corresponds to
carrier gas distribution, see Fig. 8(b). The observed differences between
studied samples could be explained in the carrier gas density and ion-carriers
interactions frame. The magnetic ions incorporated in the center of the UW1500
sample have a higher overlap with carrier gas wavefunctions than the ions on
the sides of UW1499 or UW1498. Possibly, that makes the sample UW1500 more
sensitive to the effects corresponding to the free carriers than the other
samples. Moreover, the asymmetry between samples UW1499 and UW1498 discussed
in the previous part can also affect the SLR time in the way we observe it.
Namely, the relaxation for the UW1498 is slower – as the overlap of the holes
and ions for this sample is lower due to asymmetry caused by the electric
field.

##  V Conclusions

Our work shows that the ODMR technique is a powerful tool for studying the
distribution of the strain in the nanostructures along the growth axis. It
also reveals interactions between magnetic ions and charge carriers. Both
strain and interactions with carriers can affect SLR, which may be essential
for future applications. With the ODMR technique, we performed measurements on
the series of samples with Mn2+ ions incorporated in the different regions of
the QWs – on the sides – near the barrier layers, and in the center. In
addition, we found that we can probe the different layers of the well along
the growth axis using the excited excitonic states.

We do not see any indications that the deformation of the crystal lattice is
changing along the growth axis. Therefore, we conclude that the growth of the
QW layers is homogeneous.

Additionally, the carriers play an essential role in the presented samples.
Even if the hole density in the QWs seems to be low, the introduced electric
field affects the shape of the confined wavefunctions and results in a
consecutive change of their overlap with magnetic ions. Effects related to
free charges are more pronounced for low magnetic fields – which is in
agreement with the previous studies Scherbakov  _et al._ (2001a); König  _et
al._ (2000). We conclude that the overlap between magnetic ions and hole gas
is the main factor varying across this series of QWs. Of special interest is
the case of assymetrically doped QWs, which evidence the importance of the
electric field across the structure. Our observations prove the potential to
cotrol the SLR time e.g., by deterministic tuning of the electric field in
gated samples.

##  VI Acknowledgements

This work was supported by the Polish National Science Centre under decisions
DEC-2016/23/B/ST3/03437, DEC-2015/18/E/ST3/00559, DEC-2020/38/E/ST3/00364,
DEC-2020/39/B/ST3/03251.

## References

  * Deveaud  _et al._ (1984) B. Deveaud, J. Emery, A. Chomette, B. Lambert, and M. Baudet, Applied Physics Letters 45, 1078 (1984).
  * Gaj  _et al._ (1994) J. A. Gaj, W. Grieshaber, C. Bodin-Deshayes, J. Cibert, G. Feuillet, Y. Merle d’Aubigné, and A. Wasiela, Physical Review B 50, 5512 (1994).
  * Regreny  _et al._ (1987) A. Regreny, P. Auvray, A. Chomette, B. Deveaud, G. Dupas, J. Emery, and A. Poudoulec, Revue de Physique Appliquée 22, 273 (1987).
  * Grieshaber  _et al._ (1994) W. Grieshaber, C. Bodin, J. Cibert, J. Gaj, Y. Merle d’Aubigné, A. Wasiela, and G. Feuillet, Applied Physics Letters 65, 1287 (1994).
  * Kossacki  _et al._ (1995) P. Kossacki, Nguyen The Khoi, J. A. Gaj, G. Karczewski, T. Wojtowicz, J. Kossut, and K. V. Rao, Solid State Communications 94, 439 (1995).
  * Scherbakov  _et al._ (2001a) A. V. Scherbakov, D. R. Yakovlev, A. V. Akimov, I. A. Merkulov, B. König, W. Ossau, L. W. Molenkamp, T. Wojtowicz, G. Karczewski, G. Cywinski, and J. Kossut, Physical Review B 64, 155205 (2001a).
  * König  _et al._ (2000) B. König, I. A. Merkulov, D. R. Yakovlev, W. Ossau, S. M. Ryabchenko, M. Kutrowski, T. Wojtowicz, G. Karczewski, and J. Kossut, Physical Review B 61, 16870 (2000).
  * Lambe and Kikuchi (1960) J. Lambe and C. Kikuchi, Physical Review 119, 1256 (1960).
  * Wolos  _et al._ (2004) A. Wolos, M. Palczewska, M. Zajac, J. Gosk, M. Kaminska, A. Twardowski, M. Bockowski, I. Grzegory, and S. Porowski, Physical Review B 69, 115210 (2004).
  * Bogucki  _et al._ (2022) A. Bogucki, M. Goryca, A. Łopion, W. Pacuski, K. E. Połczyńska, J. Z. Domagała, M. Tokarczyk, T. Fąs, A. Golnik, and P. Kossacki, Physical Review B 105, 075412 (2022).
  * van Gisbergen  _et al._ (1993) S. J. C. H. M. van Gisbergen, M. Godlewski, R. R. Gała̧zka, T. Gregorkiewicz, C. A. J. Ammerlaan, and N. T. Khoi, Physical Review B 48, 11767 (1993).
  * Ivanov  _et al._ (2001) V. Ivanov, M. Godlewski, S. Yatsunenko, A. Khachapuridze, M. Li, and Z. Gołacki, Acta Physica Polonica A 100, 351 (2001).
  * Ivanov  _et al._ (2007) V. Y. Ivanov, M. Godlewski, D. Yakovlev, S. Ryabchenko, G. Karczewski, and A. Waag, Physica Status Solidi (A) 204, 174 (2007).
  * Tolmachev  _et al._ (2020) D. O. Tolmachev, V. Yu. Ivanov, D. R. Yakovlev, E. V. Shornikova, B. Witkowski, S. Shendre, F. Isik, S. Delikani, H. Volkan Demir, and M. Bayer, Nanoscale 12, 21932 (2020).
  * Shornikova  _et al._ (2020) E. V. Shornikova, D. R. Yakovlev, D. O. Tolmachev, V. Y. Ivanov, I. V. Kalitukha, V. F. Sapega, D. Kudlacik, Y. G. Kusrayev, A. A. Golovatenko, S. Shendre, S. Delikanli, H. V. Demir, and M. Bayer, ACS Nano 14, 9032 (2020).
  * Gaj and Kossut (2010) J. A. Gaj and J. Kossut, eds.,  _Introduction to the Physics of Diluted Magnetic Semiconductors_, Springer Series in Materials Science (Springer-Verlag, Berlin Heidelberg, 2010).
  * Maślana  _et al._ (2003) W. Maślana, P. Kossacki, M. Bertolini, H. Boukari, D. Ferrand, S. Tatarenko, J. Cibert, and J. A. Gaj, Applied Physics Letters 82, 1875 (2003).
  * Goryca and Bogucki (2021) M. Goryca and A. Bogucki, “Sample holder for measurements of optically detected magnetic resonance,” "Sample holder patent WO2021064687A1" (2021).
  * Waag  _et al._ (1993) A. Waag, H. Heinke, S. Scholl, C. R. Becker, and G. Landwehr, Journal of Crystal Growth 131, 607 (1993).
  * Gaj  _et al._ (1979) J. Gaj, R. Planel, and G. Fishman, Solid State Communications 29, 435 (1979).
  * Cibert  _et al._ (1993) J. Cibert, R. André, C. Bodin, L. S. Dang, G. Feuillet, and P.-H. Jouneau, Physica Scripta 1993, 487 (1993).
  * Vanelle  _et al._ (1996) E. Vanelle, A. Alexandrou, J.-P. Likforman, D. Block, J. Cibert, and R. Romestain, Physical Review B 53, R16172 (1996).
  * Grieshaber  _et al._ (1996) W. Grieshaber, A. Haury, J. Cibert, Y. M. d’Aubigné, A. Wasiela, and J. Gaj, Physical Review B 53, 4891 (1996).
  * Kossacki  _et al._ (1999) P. Kossacki, J. Cibert, D. Ferrand, Y. M. d’Aubigné, A. Arnoult, A. Wasiela, S. Tatarenko, and J. Gaj, Physical Review B 60, 16018 (1999).
  * Kossacki (2003) P. Kossacki, Journal of Physics: Condensed Matter 15, R471 (2003).
  * Kossacki  _et al._ (2004) P. Kossacki, H. Boukari, M. Bertolini, D. Ferrand, J. Cibert, S. Tatarenko, J. A. Gaj, B. Deveaud, V. Ciulin, and M. Potemski, Physical Review B 70 (2004), 10.1103/PhysRevB.70.195337.
  * Łopion  _et al._ (2020) A. Łopion, A. Bogucki, K. E. Połczyńska, W. Pacuski, A. Golnik, T. Kazimierczuk, and P. Kossacki, Journal of Electronic Materials (2020), 10.1007/s11664-020-08181-z.
  * Qazzaz  _et al._ (1995) M. Qazzaz, G. Yang, S. H. Xin, L. Montes, H. Luo, and J. K. Furdyna, Solid State Communications 96, 405 (1995).
  * Strutz  _et al._ (1992) T. Strutz, A. M. Witowski, and P. Wyder, Physical Review Letters 68, 3912 (1992).
  * Goryca  _et al._ (2015) M. Goryca, M. Koperski, T. Smoleński, Ł. Cywiński, P. Wojnar, P. Plochocka, M. Potemski, and P. Kossacki, Physical Review B 92 (2015), 10.1103/PhysRevB.92.045412.
  * Strutz  _et al._ (1993) T. Strutz, A. Witowski, and P. Wyder, Review of Scientific Instruments 64, 1853 (1993).
  * Akimov  _et al._ (2006) A. Akimov, A. Scherbakov, D. Yakovlev, I. Merkulov, M. Bayer, A. Waag, and L. Molenkamp, Physical Review B 73, 165328 (2006).
  * Scherbakov  _et al._ (2001b) A. Scherbakov, A. Akimov, D. Yakovlev, W. Ossau, L. Molenkamp, S. Tatarenko, and J. Cibert, Solid State Communications 120, 17 (2001b).


# A free from local minima algorithm for training regressive MLP neural
networks

Augusto MONTISCI<EMAIL_ADDRESS>
Electrical and Electronic Engineering Dept.
University of Cagliari
Italy

###### Abstract

In this article an innovative method for training regressive MLP networks is
presented, which is not subject to local minima. The Error-Back-Propagation
algorithm, proposed by William-Hinton-Rummelhart, has had the merit of
favouring the development of machine learning techniques, which has permeated
every branch of research and technology since the mid-1980s. This
extraordinary success is largely due to the black-box approach, but this same
factor was also seen as a limitation, as soon more challenging problems were
approached. One of the most critical aspects of the training algorithms was
that of local minima of the loss function, typically the mean squared error of
the output on the training set. In fact, as the most popular training
algorithms are driven by the derivatives of the loss function, there is no
possibility to evaluate if a reached minimum is local or global. The algorithm
presented in this paper avoids the problem of local minima, as the training is
based on the properties of the distribution of the training set, or better on
its image internal to the neural network. The performance of the algorithm is
shown for a well-known benchmark.

Keywords: Multi Layer Perceptrons, Training algorithm, Local minima, Internal
representation, Non-aggregate loss function

##  1 Introduction

Even if Machine Learning (ML) includes a multiciplity of paradigms, much
different among them, most part can be considered as an evolution of Error
Backpropagation Algorithm (EBP) of Multi Layer Perceptron (Rumelhart et al.
(1986)). The merit of this algorithm consists in the fact that for the first
time it was possible to train networks with an intermediate layer, and
therefore reproduce non-linear input-output relationships. This was as true
for classification problems as it was for regression problems. Concerning the
latter ones, Cybenko (1989) had demonstrated that a Perceptron with a single
hidden layer is a universal approximator, however leaving open the problem of
determining both the number of neurons needed to solve a specific problem, and
how to determine the connection weights. The EBP algorithm offered a tool to
determine the value of the parameters, while the determination of the optimal
number of neurons is still an open problem. Successively, methods have been
presented to address both questions (Delogu et al. (2008), Ploj (2014), Ploj
et al. (2014), Fernández-Delgado et al. (2014), Fernández-Delgado et al.
(2011), Ploj et al. (2011), Curteanu and Cartwright (2011), Carcangiu et al.
(2009a)), but the EBP paradigm, with an important set of variations, still
represents the standard of machine learning. This paradigm consists in finding
the minimum of a loss function, which is typically given by the output mean
squared error with respect to the target value. All the minimization
techniques also developed in contexts other than that of ML have been proposed
to solve this problem, but the standard is represented by the use of first and
second order descent methods (Marquardt (1963)), in whose category the EBP
itself falls. First order algorithms, such as EBP, have made a comeback with
the advent of Deep Learning (LeCun et al. (2015)), as the huge number of
parameters makes second order methods impractical, even in cases in which
approximate expressions of the Hessian are used. The methods based on the
derivatives of the loss function have the advantage of being simple to
implement, but lack a criterion that allows to establish whether a
stationarity point of the loss function represents a global minimum or not.
However, there is another limitation that derives from this type of black-box
approach. In fact, the fact of considering the loss function as the only
indicator of the network performance leads to not distinguishing the different
functions of the single parts of the network, as well as the relevance of
single examples of the training set.

Figure 1: Multilayer Perceptron scheme

As a premise to the description of the algorithm that will be presented in the
next sections, it is worth rethinking the structure of the MLP network, going
beyond the black-box interpretation. First of all, the need for the hidden
layer derives from the fact that the input-output relationship is nonlinear.
If this were not the case, a linear hidden layer would be completely useless,
because the cascade of the two connecting layers would be completely
equivalent to their product. Furthermore, since the algebraic relationship
between the output of the hidden layer, called Feature Space (Cortes and
Vapnik (1995)) and the output of the neural network is linear, we can think
that the function of the hidden layer is to make linear the relationship with
the output. In other words, for each output neuron, the points of the training
set in the product space $\\{$feature space$\\}\times\\{$output space$\\}$
must be coplanar. In fact, any residual nonlinearity downstream of the feature
space could not be corrected by the last layer of connections. Hence also the
fact that there is no reason to include in the loss function the weights of
the connections preceding the output, since given the points in the feature
space, the optimal set of weights is the one that corresponds to the linear
regression plane. Another important consideration concerns the different role
of the hidden layer and the output layer. Indeed, if the former represents the
degrees of freedom of the network, because as the size of the feature space
increases, it becomes easier to make the points of the training set coplanar
in the product space $\\{$feature space$\\}\times\\{$output space$\\}$, the
latter constitutes a constraint, because the same distribution in the feature
space must satisfy the coplanarity constraint with respect to different
outputs.

This paper presents a training algorithm based on this interpretation of the
algebraic structure of the MLP. For sake of simplicity and without prejudice
to the general validity of the results, in this article the case of regressive
MLP networks with only one output (MISO) will be treated. The paper is
structured as follows: in the first section, the analytical basis of the
procedure is presented. The second section briefly describes the choices
adopted for implementation in the Matlab environment. The third shows the
results obtained with a well-known benchmark. At the end, conclusions are
provided.

##  2 Description of the algorithm

$\begin{cases}\mathbf{W\cdot x+d=y}\\\ \mathbf{h=\sigma(y)}\\\ \mathbf{V\cdot
h+b=u}\end{cases}$ (1)

In this section the training algorithm will be described, referring to the
scheme shown in Fig. 1. The choice of having only one output serves to
simplify the discussion, but does not affect the generality of the problem.
The MLP network implements an algebraic structure represented by the system of
equations (1), where $\boldsymbol{x}$ is the input to the network,
$\boldsymbol{W}$ is the weight matrix of the first connection layer,
$\boldsymbol{d}$ is the bias of the first layer, $\boldsymbol{y}$ is the input
vector to the hidden layer, $\sigma$ is the activation function of the hidden
layer, $\boldsymbol{h}$ is the image of the input in the feature space,
$\boldsymbol{V}$ is the vector of the weights of the connections with the
output, $\boldsymbol{b}$ is the bias of the second layer. The first layer of
connections linearly transforms the points from the input space of the network
into the input space of the hidden layer. The hidden layer activation function
is the only nonlinearity of the network. Through the first two transformations
of the (1), the distribution of the training set points in the input space is
transformed into a distribution in the feature space. The first layer of
connections must ensure that the points of the training set in the product
space $\\{$feature space$\\}\times\\{$output space$\\}$ are all coplanar. If
this does not happen, the second layer of connections will not be able to
compensate for these errors. The solution that minimizes least squares is the
regression plan. For this reason, even if a conventional training algorithm is
applied, it is not convenient to include the weights of the second layer of
connections within the loss function. In fact, it is better to define the loss
as a function of the weights of the first layer and, by freezing this, to
define the second layer by calculating the regression plane.

The algorithm presented in this paper, instead of finalizing the search to
minimize a loss function, aims to construct the feature space in such a way as
to ensure the coplanarity of the points, or equivalently that the level curves
of the function in the input space are transformed into straight lines in the
feature space (Fig. 2).

Figure 2: Linearizzation performed by the hidden layer

The condition for linearization can be expressed analytically as follows:

$\mathbf{s^{T}\cdot\sigma(W\cdot
x_{i}+d)}=\textit{a}\cdot\mathbf{u_{i}}\quad\forall\textit{i}=1\dots\textit{N}$
(2)

where $\boldsymbol{s}$ is a vector of linear combination coefficients of the
points in the feature space, while $a$ is a scalar of arbitrary value, which
determines the distance between the straight lines of level, and consequently
the slope of the regression plan. Finally, $N$ is the number of points in the
training set. Dividing the equation (2) by $a$ eliminates the uncertainty of
the solution:

$\mathbf{\hat{s}^{T}\cdot\sigma(W\cdot
x_{i}+d)}=\mathbf{u_{i}}\quad\forall\textit{i}=1\dots\textit{N}$ (3)

where $\hat{\boldsymbol{s}}$ is the normalized variable. Globally the equation
(3) is a system of $N$ nonlinear equations in the variables $\boldsymbol{W}$,
$\boldsymbol{d}$, $\boldsymbol{\hat{s}}$. Fixing the values of
$\boldsymbol{W}$ and $\boldsymbol{d}$ the same system becomes a system of
linear equations overdetermined in the single variable $\hat{\boldsymbol{s}}$.

The proposed algorithm starts from a first order approximation of the first
member of the equation (3). Let $\boldsymbol{W}_{0}$, $\boldsymbol{d}_{0}$,
$\boldsymbol{\hat{s}}_{0}$ be an initial solution of the system (3), and
$\delta\boldsymbol{W}$, $\delta\boldsymbol{d}$, $\delta\boldsymbol{\hat{s}}$
the respective increments of the three variables that provide the solution:

$\mathbf{(\hat{s}_{0}+\delta s)^{T}\cdot\sigma[(W_{0}+\delta W)\cdot
x_{i}+(d_{0}+\delta
d)]}=\mathbf{u_{i}}\quad\forall\textit{i}=1\dots\textit{N}$ (4)

Now replace the left-hand side of the equation with a first-order
approximation of the incremented function:

$\mathbf{(\hat{s}_{0}+\delta
s)^{T}\cdot[\sigma_{0}+\nabla\sigma(W,d)\odot\delta_{W,d}\cdot
1]}\approx\mathbf{u_{i}}\quad\forall\textit{i}=1\dots\textit{N}$ (5)

where $\boldsymbol{\sigma_{0}}$ is the vector in the feature space
corresponding to a zero increment, $\nabla\sigma(\boldsymbol{W,d})$ is the
gradient calculated with respect to the weights of the first layer
connections, $\delta_{\boldsymbol{W,d}}$ is the increment of the variables,
$\odot$ is the elementwise multiplication, $\boldsymbol{1}$ is a column vector
of only 1s. Reordering the equation and neglecting the infinitesimals of
higher order we get:

$\mathbf{\sigma_{0}^{T}\cdot\delta
s+\hat{s}_{0}^{T}\cdot\nabla\sigma(W,d)\odot\delta_{W,d}\cdot
1}\approx\mathbf{u_{i}-\hat{s}_{0}^{T}\cdot\sigma_{0}}\quad\forall\textit{i}=1\dots\textit{N}$
(6)

The equation (6) defines a linear, overdetermined system of equations, whose
unknowns are the increments of the parameters. By combining (3) with (6) it is
possible to implement an iterative procedure which provides the weights of the
first layer of connections. Taking into account the fact, as previously
indicated, that the second layer of connections is uniquely determined once
the weights of the first layer have been assigned, the iterative procedure
actually completes the training procedure.

The procedure takes place as follows. Given the initial value of the weights
of the first layer of connections, the system (3) is resolved into the single
variable $\boldsymbol{\hat{s}}$. Since the system is overdetermined, it can
only be solved by minimizing the root mean square error (Jódar et al. (1991)).
Using the obtained solution, together with the same weights of the first layer
used previously, the coefficients and the known term of the system (6) are
obtained. This system is also overdetermined, and therefore it too will have
to be solved according to least squares. Only the increments of the weights of
the first layer $\delta_{\boldsymbol{W,d}}$ of the solution of the system are
used, while $\boldsymbol{\hat{s}}$ is updated by solving the system (3) after
updating the coefficients based on the new values of $\mathbf{W}$ and
$\mathbf{d}$. The iterative process ends when the error in the system (3) is
less than a pre-set threshold value. Alternatively, since the second layer of
weights are uniquely determined by the first layer of weights, the error value
of the current solution can be computed directly on the neural network output.
Since a first-order approximation is being used, the question remains open as
to what is the appropriate step length to take at each iteration. For this
reason, the Line Search method is adopted (Bazaraa et al. (2013)), sampling
the segment from the current solution to the increased one, and evaluating for
each point the value of the output error. At each iteration, the assigned
increment will be the one that corresponds to the minimum error.

As can be deduced, the number of operations to be performed for each iteration
is much higher than that of a conventional method, but the proposed method has
the advantage of not being subject to local minima problems, since the
objective is not more that of minimizing a loss function, rather than solving
a system of equations, and therefore the loss function that drives the search
for the solution coincides with the goal of the training, rather than, as
usually happens, being constituted by an aggregate target function.

It is worth noting that the proposed method, precisely because it is not based
on the minimization of an aggregate loss function, allows the use of different
criteria in the application of the line search, such as for example the
minimization of the maximum error in the training set. The complete analysis
of the potential of the method goes beyond the objectives of this paper and
will be the subject of future developments.

##  3 Results

A well-known benchmark for testing optimization methods (Schwefel (1981)) is
used as an application example. In cases where the search for the optimal
solution is carried out with the aid of a neural network (Carcangiu et al.
(2009b)), the accuracy of the network and its generalization capability assume
fundamental importance. The general expression of the function is:

$f(\textbf{x})=418.9829\cdot d-\sum_{i=1}^{d}x_{i}\sin(\sqrt{\lvert
x_{i}\rvert})$ (7)

where $d$ is the size of the space where the function is defined. In this case
a dimension $d=3$ was used, and a range of variables [-500,500]. Fig. 3 shows
the Schwefel function used as a test. The value of the function is represented
by the color gradient. It is apparent that the function has a large number of
maxima and minima, which implies that a large number of hidden neurons will be
required to adequately fit it.

Figure 3: Schwefel’s function defined in d=3

An MLP network with 200 hidden nodes has been trained for 4000 epochs. Fig. 4
shows the trend of the mean squared error during training. Note that the error
in the diagram refers to the output data previously normalized between -1 and
1. The total number of training examples was equal to 5152. As can be seen, in
the first iterations there is a considerable reduction of the error, while
subsequently the descent speed is considerably reduced both in absolute and
relative terms. In some segments it can be noted that the average error
appears constant, and in some cases it can be slightly increasing, but this
does not prevent, after a suitable number of iterations, from resuming the
decreasing trend. This means that although the mean (squared) error is
constant, in reality the iterative process is modifying the configuration of
the weights, which allows the network to identify regions in the parameter
space where it is possible to obtain performance improvements. As can be seen,
the error value after 4000 epochs had not yet stabilized, although the
decrease had become much slower.

As mentioned earlier, it is possible to redefine the search criterion based on
the maximum error in the training set rather than the mean value, but this,
like other issues affecting convergence acceleration, is beyond the scope of
this paper and will be the subject of future works.

Figure 4: Performance diagram of the training

##  4 Some implementation notes

The algorithm was developed in Matlab environment. This has led, as far as
possible, to exploit matrix calculation, rather than resorting to for loops.
This applies in particular to the construction, at each iteration, of the
matrix of the system coefficients, which has made it possible to enormously
contain the time required for a single iteration.

It is also worth making some considerations on the application of the Line
Search, with which the length of the step to be taken at each iteration is
established. The number of samples that are explored for each iteration is a
macro-parameter that must be set by the designer, and significantly affects
the calculation time of the single iteration. In the case under examination,
the increment obtained as a solution of the linear system is divided into 1000
points distributed with a logarithmic law, so as to better detail the trend of
the function for small increments. It may happen that the minimum of the error
corresponds to the first sample, which implies that the minimum could
correspond to a smaller increment than the minimum foreseen. In this case the
search for the minimum is further refined using the bisection method applied
to the interval between 0 and the smallest increment. If the minimum falls in
correspondence with a zero increase, the training procedure stops.

The procedure requires solving a large number of linear systems. The use of
the Matlab library functions in the case where the matrix of the coefficients
is ill-conditioned involves difficulties, both in terms of calculation times
and instability of the iterative process. For this reason, also taking into
account that the system is approximated in its definition, we use a solution
method based on successive projections (Cannas et al. (2012)), rather than on
the inversion of the coefficient matrix. The obtained solution is anyway
affected by the misconditioning of the coefficient matrix, and this generally
results in a less precise solution, but the algorithm has the merit of quickly
providing a reasonable solution.

Again regarding the problem of the misconditioning of coefficient matrices, it
has been observed that a large number of examples of the training set helps to
mitigate the problem. Increasing the training set does not significantly
affect the calculation time, if not indirectly due to memory occupation. It
has also been observed that the calculation time of the single iteration does
not significantly depend on the number of hidden neurons, therefore it is
observed that as the number of hidden neurons increases, the total training
time decreases, because the number of epochs necessary to reach the target
error is smaller.

As stopping criterion a double check was adopted on the number of epochs and
on the minimum variation of the error between two successive epochs. In all
the tested cases, irrespective of the number of examples, the algorithm
stopped when the maximum number of epochs was reached. In future developments
it is planned to use a stopping criterion based on the stabilization of the
weights of the connections, rather than on the outputs. This will prevent the
procedure from stopping due to a momentary stabilization of the error trend,
like the one found in the example in the previous section. If, in fact, it is
observed that the weights of the connections continue to change, even if this
does not translate into a reduction in the average error of the network, the
training must still proceed. Analyzing weight trends during training offers
the possibility of accelerating learning. If, in fact, this trend is regular,
it is generally possible to extrapolate the values to predict which value it
is tending to. This possibility will be the subject of future developments.

##  5 Conclusion

This paper describes a training algorithm for Multi Layer Perceptron neural
networks whose main feature is that it is not subject to the problem of local
minima, typical of conventional algorithms, in which the search for the
minimum of the loss function is driven by derivatives of the first and second
order of the loss function. Even at a preliminary implementation, the
algorithm shows good performance and limited computational complexity as both
the training set and the network size increase, which makes it a possible
candidate for the treatment of big data problems using neural networks. The
performance of the algorithm has been tested on a benchmark that is difficult
to treat with the most common training algorithms, obtaining very encouraging
results. Finally, the paper introduces lines of development that will be the
subject of future publications.

## References

  * Bazaraa et al. (2013) M.S. Bazaraa, H.D. Sherali, and C.M. Shetty.  _Nonlinear programming: theory and algorithms_.  John wiley & sons, 2013.
  * Cannas et al. (2012) B. Cannas, S. Carcangiu, A. Fanni, R. Forcinetti, A. Montisci, G. Concu, et al.  Hops: A new tomographic reconstruction algorithm for non destructive acoustic testing of concrete structures.  _Advances in civil engineering and building materials_ , 831:259–263, 2012.
  * Carcangiu et al. (2009a) S. Carcangiu, A. Fanni, and A. Montisci.  A constructive algorithm of neural approximation models for optimization problems.  _COMPEL-The international journal for computation and mathematics in electrical and electronic engineering_ , 28(5):1276–1289, 2009a.
  * Carcangiu et al. (2009b) S. Carcangiu, A. Fanni, and A. Montisci.  _Multi objective optimization algorithm based on neural networks inversion_ , pages 744–751.  2009b.
  * Cortes and Vapnik (1995) C. Cortes and V. Vapnik.  Support-vector networks.  _Machine Learning_ , 20(3):273 – 297, 1995.
  * Curteanu and Cartwright (2011) S. Curteanu and H. Cartwright.  Neural networks applied in chemistry. i. determination of the optimal topology of multilayer perceptron neural networks.  _Journal of Chemometrics_ , 25(10):527 – 549, 2011.
  * Cybenko (1989) G. Cybenko.  Approximation by superpositions of a sigmoidal function.  _Mathematics of Control, Signals, and Systems_ , 2(4):303 – 314, 1989.
  * Delogu et al. (2008) R. Delogu, A. Fanni, and A. Montisci.  Geometrical synthesis of mlp neural networks.  _Neurocomputing_ , 71(4-6):919 – 930, 2008.
  * Fernández-Delgado et al. (2011) M. Fernández-Delgado, J. Ribeiro, E. Cernadas, and S.B. Ameneiro.  Direct parallel perceptrons (dpps): Fast analytical calculation of the parallel perceptrons weights with margin control for classification tasks.  _IEEE Transactions on Neural Networks_ , 22(11):1837 – 1848, 2011.
  * Fernández-Delgado et al. (2014) M. Fernández-Delgado, E. Cernadas, S. Barro, J. Ribeiro, and J. Neves.  Direct kernel perceptron (dkp): Ultra-fast kernel elm-based classification with non-iterative closed-form weight calculation.  _Neural Networks_ , 50:60 – 71, 2014.
  * Jódar et al. (1991) L. Jódar, A.G. Law, A. Rezazadeh, and J.H. Weston.  Computations for the moore penrose and other generalized inverses.  _Congressus Numerantium_ , pages 57–57, 1991.
  * LeCun et al. (2015) Y. LeCun, Y. Bengio, and G. Hinton.  Deep learning.  _nature_ , 521(7553):436–444, 2015.
  * Marquardt (1963) D. W. Marquardt.  An algorithm for least-squares estimation of nonlinear parameters.  _Journal of the society for Industrial and Applied Mathematics_ , 11(2):431–441, 1963.
  * Ploj (2014) B. Ploj.  Optimization for multi layer perceptron: Without the gradient, 2014.
  * Ploj et al. (2011) B. Ploj, M. Zorman, and P. Kokol.  Border pairs method - constructive mlp learning classification algorithm.  _Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics)_ , 6943 LNAI:297 – 307, 2011.
  * Ploj et al. (2014) B. Ploj, R. Harb, and M. Zorman.  Border pairs method-constructive mlp learning classification algorithm.  _Neurocomputing_ , 126:180 – 187, 2014.
  * Rumelhart et al. (1986) D.E. Rumelhart, G. E. Hinton, and R. J. Williams.  Learning representations by back-propagating errors.  _Nature_ , 323(6088):533 – 536, 1986.
  * Schwefel (1981) H.-P. Schwefel.  _Numerical optimization of computer models_.  John Wiley & Sons, Inc., 1981.


# Factorials and Legendre’s three-square theorem: II

Rob Burns

###### Abstract

Let $\bar{S}$ denote the set of integers $n$ such that $n!$ cannot be written
as a sum of three squares. Let $\bar{S}(n)$ denote $\bar{S}\cap[1,n]$. We
establish an exact formula for $\bar{S}(2^{k})$ and show that
$\bar{S}(n)=1/8*n+\mathcal{O}(\sqrt{n})$. We also list the lengths of gaps
appearing in $\bar{S}$. We make use of the software package Walnut to
establish these results.

##  1 Introduction

In a previous paper [1], we constructed an automaton accepting integers $n$
such that $n!$ cannot be written as a sum of three squares. In the current
paper we use this automaton and the software package Walnut to explore the
sets

$S:=\\{n:n!\,\,\text{is a sum of three squares}\\},$
$\bar{S}:=\\{n:n!\,\,\text{is not a sum of three squares}\\},$ $S(n):=\\{x\leq
n:x!\,\,\text{is a sum of three squares}\\}.$

and

$\bar{S}(n):=\\{x\leq n:x!\,\,\text{is not a sum of three squares}\\}.$

In particular, we establish a formula for $S(2^{k})$, or equivalently,
$\bar{S}(2^{k})$, and show that $\bar{S}(n)=1/8*n+\mathcal{O}(\sqrt{n})$. We
also list the lengths of possible gaps in $S$ and its complement $\bar{S}$.

We mention here some similar previous results. In [2], Deshouillers and Luca
established the two equivalent bounds

$S(n)=7n/8+\mathcal{O}(n^{2/3})\,\,\text{ and
}\,\,\bar{S}(n)=x/8+\mathcal{O}(n^{2/3}).$

This bound was improved by Hajdu and Papp [7] to

$\bar{S}(n)=x/8+\mathcal{O}(n^{1/2}\log^{2}n).$

Hajdu and Papp [7] proved that the largest gap between consecutive elements of
$\bar{S}$ is $42$. In theorem 3.1, we list the lengths of all gaps that occur
between consecutive elements of $S$ and $\bar{S}$.

We make use of the software package Walnut. Hamoon Mousavi, who wrote the
program, has provided an introductory article [8]. Papers that have used
Walnut include [9], [12], [10], [5], [3]. Further resources related to Walnut
can be found at Jeffrey Shallit’s page

    https://www.cs.uwaterloo.ca/~shallit/papers.html   .

The free open-source mathematics software system SageMath [14] was essential
for the matrix algebra and factoring of polynomials.

In the next section, we describe the automaton that will be used to examine
$\bar{S}$ and $\bar{S}(x)$.

##  2 The automaton

We first summarise the conditions satisfied by integers in $\bar{S}$. These
come from [1] but were previously established by Granville and Zhu in 1990
[6]. Let $n\in\mathbb{N}$, with binary representation given by $n=\sum_{k\geq
0}a_{k}2^{k}$, where all but finitely many $a_{i}$ are zero. Let $\gamma$ be
the highest power of $2$ dividing $n!$. For $i\in\\{3,5\\}$, define
$\alpha_{i}=\alpha_{i}(n)$ by

$\displaystyle\alpha_{3}=\alpha_{3}(n)$ $\displaystyle:=\\#\Big{\\{}k\geq
0:\sum_{i=k}^{k+2}a_{i}2^{i-k}\in\\{3,4\\}\Big{\\}}$ (1)
$\displaystyle\alpha_{5}=\alpha_{5}(n)$ $\displaystyle:=\\#\Big{\\{}k\geq
0:\sum_{i=k}^{k+2}a_{i}2^{i-k}\in\\{5,6\\}\Big{\\}}.$ (2)

######  Theorem 2.1.

Let $n\in\mathbb{N}$ have binary representation $n=\sum_{k\geq 0}a_{k}2^{k}$,
where all but finitely many $a_{i}$ are zero. If $\gamma$ is the highest power
of $2$ dividing $n!$, then $n!=2^{\gamma}Z$, where $Z$ satisfies

$\displaystyle Z\equiv 3^{\alpha_{3}(n)}(-1)^{\alpha_{5}(n)}\pmod{8}$

######  Corollary 2.2.

If $n\in\mathbb{N}$, then $n!$ cannot be written as a sum of three squares if
and only if $\gamma$ and $\alpha_{3}$ are even and $\alpha_{5}$ is odd.

According to the above two results, whether $n!$ can be written as a sum of
three squares is determined by the triple
$(\gamma\pmod{2},\alpha_{3}\pmod{2},\alpha_{5}\pmod{2})$. For convenience, we
will give this triple a name. For integers $n$, define $\Theta(n)$ by

$\Theta(n):=(\,\,\gamma\pmod{2},\,\,\alpha_{3}\pmod{2},\,\,\alpha_{5}\pmod{2}\,\,).$

The components of $\Theta(n)$ are

$\Theta_{1}(n):=\gamma\pmod{2},\,\,\,\Theta_{2}(n):=\alpha_{3}\pmod{2},\,\,\,\Theta_{3}(n):=\alpha_{5}\pmod{2}.$
(3)

Then, $n\in\bar{S}$ if and only if $\Theta(n)=(0,0,1)$.

We present here three automata which take as input the binary digits of an
integer $n$, starting with the least significant digit. The automata were
produced by Walnut based on the instructions provided in [1]. The first
automaton, which we call gamma, is pictured in figure 1. It accepts integers
$n$ if and only if $\gamma(n)$ is even. The second automaton, which we call
a3, is pictured in figure 2. It accepts integers $n$ if and only if
$\alpha_{3}(n)$ is even. The third automaton, which we call a5, is pictured in
figure 3. It accepts integers $n$ if and only if the $\alpha 5(n)$ is even.

Figure 1: Automaton for the parity of $\gamma$ Figure 2: Automaton for the
parity of $\alpha_{3}$ Figure 3: Automaton for the parity of $\alpha_{5}$

For any triple of binary digits $(x,y,z)$, we can use Walnut to create an
automaton which accepts $n$ if and only if $\Theta(n)=(x,y,z)$. We are
particularly interested in determining when $\Theta(n)=(0,0,1)$. We create
this automaton using the Walnut command:

    def factauto "?lsd_2 $gamma(n) & $a3(n) & ~$a5(n)":

The resulting automaton has 33 states and is pictured in figure 4.

Figure 4: Automaton accepting $n$ if $n!$ cannot be written as a sum of $3$
squares

##  3 Gaps and runs in $S$ and $\bar{S}$

In this section we investigate the gaps and runs in the sets $S$ and
$\bar{S}$. A gap $g$ exists in $S$ (or $\bar{S}$) if there are integers $x$
and $y$ which are consecutive elements of $S$ and such that $y-x=g$. A run of
length $r$ exists in $S$ (or $\bar{S}$) if $S$ contains $x,x+1,\dots,x+r-1$
but $x-1$ and $x+r$ are not in $S$. Since $S\cup\bar{S}=\mathbb{N}$, gaps in
$S$ correspond to runs in $\bar{S}$ and vice versa. For example, if $S$
contains a run of size $r$, then $\bar{S}$ contains a gap of length $r+1$.

In [7], Hajdu and Papp proved that the largest gap between consecutive
elements of $\bar{S}$ is $42$.

######  Theorem 3.1.

The gap sizes in $S$ are $\\{1,2,3,4\\}$. The gap sizes in $\bar{S}$ are

$\\{1,2,3,\dots,23,25,26,28,30,31,33,34,35,37,38,42\\}.$

###### Proof.

We firstly deal with $\bar{S}$. Recall the automaton "factauto" which was
defined previously. It accepts integers $n$ such that $n!$ cannot be written
as a sum of three squares. We use Walnut to define an automaton, which we call
"gaps":

⬇

def gaps "?lsd_2 $factauto(n) & $factauto(n+r) & (Aj (j < r-1) =>
~$factauto(n+ j + 1))":

This automaton takes as input the integer pair $(n,r)$ and reaches an
accepting state if and only if $n\in\bar{S}$ and the next integer in $\bar{S}$
is $n+r$, i.e. there a gap of length $r$ in $\bar{S}$ starting at $n$. Since
the automaton has 319 states, we will not display it here. The following
Walnut command creates an automaton which accepts an integer $r$ if and only
if there is a gap of length $r$ somewhere in $\bar{S}$:

    eval tmp "?lsd_2 E n $gaps(n, r)":

This automaton is displayed in figure 5. A patient reader will see that the
automaton accepts the integers listed in the theorem.

We can use the same approach to find all gaps in the set $S$. We define the
automaton "sgaps" which takes as input the integer pair $(n,r)$ and reaches an
accepting state if and only if $n\in S$ and the next integer in $S$ is $n+r$,
i.e. there a gap of length $r$ in $S$ starting at $n$.

⬇

def sgaps "?lsd_2 ~$factauto(n) & ~$factauto(n+r) & (Aj (j < r-1) =>
$factauto(n+ j + 1))":

The resulting automaton has 203 states. The following Walnut command creates
an automaton which accepts an integer $r$ if and only if there is a gap of
length $r$ somewhere in $S$:

    eval tmp "?lsd_2 E n $sgaps(n, r)":

This automaton is displayed in figure 6. It is clear that it only accepts the
integers $1,2,3,4$. ∎

######  Remarks.

A gap of 42 in $\bar{S}$ commences at $n=23268$ as mentioned by Hajdu and
Papp. This is the second last gap to appear in the sequence $\bar{S}$. The
last gap to appear in $\bar{S}$ is $33$, which commences at
$n=153828=(100101100011100100)_{2}$.

Theorem 3.1 is about gaps in the set $\bar{S}$ which contains integers $n$
satisfying $\Theta(n)=(0,0,1)$. The same method can be used for other sets
defined in terms of $\Theta$ (see (9) below). Hajdu and Papp showed that the
maximum gap length in any of these sets is $42$.

Figure 5: Automaton accepting gap lengths in $\bar{S}$ Figure 6: Automaton
accepting gap lengths in $S$

##  4 A formula for $\bar{S}(2^{k})$

In this section we establish a formula for $\bar{S}(2^{k})$, The method we use
is taken from the papers [4], [11], [13]. Firstly, we use Walnut to obtain a
linear representation of $\bar{S}(n)$. Walnut produces matrices $M_{0}$ and
$M_{1}$ and vectors $v$ and $w$ such that, if $n$ has binary representation
$(n_{0},n_{1},\dots,n_{k})$ where $n_{0}$ is the least significant digit, then

$\bar{S}(n)=v*M_{n_{0}}*M_{n_{1}}*\dots*M_{n_{k}}*w$ (4)

In the particular case that $n=2^{k}$, we have
$\bar{S}(2^{k})=v*M_{0}^{k}*M_{1}*w$. The theory of linear recurrences says
that there are constants $\\{c_{i}\\}$ such that $\bar{S}(2^{k})$ can be
written as

$\bar{S}(2^{k})=\sum_{i}c_{i}\lambda_{i}^{k}$ (5)

where $\\{\lambda_{i}\\}$ are the roots of the minimal polynomial of $M_{0}$.
The constants $\\{c_{i}\\}$ can be determined by calculating $\bar{S}(2^{k})$
for enough values of $k$ and plugging the values into (5) to obtain a system
of linear equations that can be solved.

######  Theorem 4.1.

For $k\geq 2$,

$\bar{S}(2^{k})=\begin{cases}2^{k-3}-2^{(k-4)/2},&\text{if
}\,\,k\equiv\\{0,2,22\\}\pmod{24}\\\ 2^{k-3}-2^{(k-3)/2},&\text{if
}\,\,k\equiv\\{1,3\\}\pmod{24}\\\ 2^{k-3}-2^{(k-5)/2},&\text{if }\,\,k\equiv
11\pmod{24}\\\ 2^{k-3}+2^{(k-4)/2},&\text{if
}\,\,k\equiv\\{14,16,18\\}\pmod{24}\\\ 2^{k-3}+2^{(k-3)/2},&\text{if
}\,\,k\equiv 17\pmod{24}\\\ 2^{k-3}-3*2^{(k-5)/2},&\text{if }\,\,k\equiv
23\pmod{24}\\\ 2^{k-3}&\text{otherwise }\end{cases}$

###### Proof.

The Walnut command

    eval sumfact n "?lsd_2 (j>=1) & (j<=n) & $factauto(j)":

produces the matrices $M_{0}$ and $M_{1}$ and the vectors $v$ and $w$
mentioned above. Since the dimensions of these objects are large, we will not
display them here. The minimal polynomial of $M_{0}$ is

$\begin{split}h&(x)=\,\,x^{20}-7*x^{19}+20*x^{18}-30*x^{17}+24*x^{16}-40*x^{14}+64*x^{13}-128*x^{11}\\\
&+160*x^{10}-384*x^{8}+960*x^{7}-1280*x^{6}+896*x^{5}-256*x^{4}\\\
&=\,\,(x-2)*(x-1)*(x-i-1)*(x+i-1)*(x-i\sqrt{2})*(x+i\sqrt{2})*\\\
&*(x+\sqrt{2})*(x+(-i/2+1/2)*\sqrt{3}+i/2-1/2)*\\\
&*(x-1/2*\sqrt{2}*\sqrt{3}+1/2*\sqrt{2})*(x+1/2*\sqrt{2}*\sqrt{3}+1/2*\sqrt{2})*\\\
&*(x+(i/2+1/2)*\sqrt{3}-i/2-1/2)*(x+(-i/2-1/2)*\sqrt{3}-i/2-1/2)*\\\
&*(x-1/2*\sqrt{2}*\sqrt{3}-1/2*\sqrt{2})*(x+1/2*\sqrt{2}*\sqrt{3}-1/2*\sqrt{2})*\\\
&*(x+(i/2-1/2)*\sqrt{3}+i/2-1/2)*(x-\sqrt{2})*x^{4}\end{split}$

All roots of $h$, apart from $2$ and $1$, are of the form
$\lambda_{i}=\sqrt{2}*(\text{a 24'th root of unity})$. We then have for $0\leq
s<24$,

$\begin{split}\bar{S}(2^{24k+s})=&\sum_{j}c_{j}\lambda_{j}^{24k+s}=c_{0}*2^{24k+s}+c_{1}+\sum_{\lambda_{j}\not
1,2}c_{j}\lambda_{j}^{24k+s}\\\
=&c_{0}*2^{24k+s}+c_{1}+2^{12k}*\sum_{\lambda_{j}\not=1,2}c_{j}\lambda_{j}^{s}\end{split}$
(6)

The values of the $\\{c_{j}\\}$ can be determined by calculating
$\bar{S}(2^{k})$ for enough $k$ using (4). We create a matrix $A$ having the
same dimension as the number of constants. The $j$’th row of $A$ is

$2^{j+2},1,\lambda_{2}^{j+2},\lambda_{3}^{j+2},\dots$

Then the vector of constants is given by

$(c_{0},c_{1},\dots)=A^{-1}*(\bar{S}(2^{2}),\bar{S}(2^{3}),\dots).$

The calculation gives $c_{0}=1/8$ and $c_{1}=0$. Since we are only interested
in the values for $c_{0}$ and $c_{1}$, an alternative approach is to use (6)
and the calculated values of, say, $\bar{S}(2^{3})=0,\bar{S}(2^{27})=16773120$
and $\bar{S}(2^{51})=281474959933440$ to create three linear equations. The
equations are:

$\begin{split}2^{3}*c_{0}+c_{1}+\tau=&0\\\
2^{27}*c_{0}+c_{1}+2^{12}*\tau=&16773120\\\
2^{51}*c_{0}+c_{1}+2^{24}*\tau=&281474959933440\end{split}$

where $\tau=\sum_{\lambda_{j}\not=1,2}c_{j}\lambda_{j}^{3}$. These can be
solved to give $c_{0}=1/8$, $c_{1}=0$ and $\tau=-1$.

Plugging the values for $c_{0}$ and $c_{1}$ into (6), we see that, for each
$0\leq s<24$, the terms $\sum_{\lambda_{j}\not=1,2}c_{j}\lambda_{j}^{s}$
satisfy

$\sum_{\lambda_{j}\not=1,2}c_{j}\lambda_{j}^{s}=\bar{S}(2^{s})-2^{s-3}$

As an example of the calculation that is required for each $s$,
$\bar{S}(2^{3})=0$ so $\sum_{\lambda_{j}\not=1,2}c_{j}\lambda_{j}^{3}=-1$.
Hence, if $n\equiv 3\pmod{24}$, then $\bar{S}(2^{n})=2^{n-3}-2^{12*\lfloor
n/24\rfloor}=2^{n-3}-2^{(n-3)/2}$. Here $\lfloor x\rfloor$ is the floor
function. A similar calculation can be made for each $s$. ∎

######  Theorem 4.2.

For $k\geq 5$,

$\bar{S}(3*2^{k})=\begin{cases}3*2^{k-3}+1,&\text{if
}\,\,k\equiv\\{7,9,13,18,19\\}\pmod{24}\\\ 3*2^{k-3}+1-2^{(k-3)/2},&\text{if
}\,\,k\equiv\\{1,3,5\\}\pmod{24}\\\ 3*2^{k-3}+1-2^{(k-2)/2},&\text{if
}\,\,k\equiv\\{0,4,10\\}\pmod{24}\\\ 3*2^{k-3}+1-2^{(k-4)/2},&\text{if
}\,\,k\equiv\\{2,6,8,12\\}\pmod{24}\\\ 3*2^{k-3}+1-2^{(k-5)/2},&\text{if
}\,\,k\equiv\\{11,23\\}\pmod{24}\\\ 3*2^{k-3}+1+2^{(k-4)/2},&\text{if
}\,\,k\equiv 14\pmod{24}\\\ 3*2^{k-3}+1+2^{(k-1)/2},&\text{if }\,\,k\equiv
15\pmod{24}\\\ 3*2^{k-3}+1+3*2^{(k-4)/2},&\text{if }\,\,k\equiv 16\pmod{24}\\\
3*2^{k-3}+1+2^{(k-3)/2},&\text{if }\,\,k\equiv 17\pmod{24}\\\
3*2^{k-3}+1-3*2^{(k-4)/2},&\text{if }\,\,k\equiv\\{20,22\\}\pmod{24}\\\
3*2^{k-3}+1-2^{(k-1)/2},&\text{if }\,\,k\equiv 21\pmod{24}\end{cases}$

###### Proof.

From (4),

$\bar{S}(3*2^{k})=v*M_{0}^{k}*M_{1}*M_{1}*w$

We can therefore use the same approach as in theorem 4.1. We have,

$\bar{S}(3*2^{k})=\sum_{i}c_{i}\lambda_{i}^{k}$ (7)

where $\\{c_{j}\\}$ are constants and $\\{\lambda_{j}\\}$ are the roots of the
minimal polynomial of $M_{0}$, which we have already calculated. The constants
$c_{0}$ and $c_{1}$ can be calculated using the values of
$\bar{S}(3*2^{3})=3,\bar{S}(3*2^{27})=50327553$ and
$\bar{S}(3*2^{51})=844424913354753$ to create three linear equations. The
equations are:

$\begin{split}2^{3}*c_{0}+c_{1}+\tau=&3\\\
2^{27}*c_{0}+c_{1}+2^{12}*\tau=&50327553\\\
2^{51}*c_{0}+c_{1}+2^{24}*\tau=&844424913354753\end{split}$

where $\tau=\sum_{\lambda_{j}\not=1,2}c_{j}\lambda_{j}^{3}$. These can be
solved to give $c_{0}=3/8$, $c_{1}=0$ and $\tau=-1$. As before, we have for
$0\leq s<24$,

$\begin{split}\bar{S}(3*2^{24k+s})=&\sum_{j}c_{j}\lambda_{j}^{24k+s}=3/8*2^{24k+s}+1+\sum_{\lambda_{j}\not
1,2}c_{j}\lambda_{j}^{24k+s}\\\
=&3*2^{24k+s-3}+1+2^{12k}*\sum_{\lambda_{j}\not=1,2}c_{j}\lambda_{j}^{s}\end{split}$
(8)

The formula for $\bar{S}(3*2^{k})$ can be deduced once (7) is used to
calculate $\sum_{\lambda_{j}\not=1,2}c_{j}\lambda_{j}^{s}$ for each $0\leq
s<24$. ∎

##  5 A bound for $\bar{S}(n)$

In this section we establish a best possible bound on $\bar{S}(n)$. We borrow
an argument from the paper by Hajdu and Papp [7].

We introduce the following sets for any triple $(x,y,z)$ of binary digits:

$\begin{split}\bar{S}_{(x,y,z)}:=&\\{n:\Theta(n)=(x,y,z)\\},\\\
\bar{S}_{(x,y,z)}(n):=&\\{r\leq n:\Theta(r)=(x,y,z)\\},\\\
\bar{S}_{(x,y,z)}(m,n):=&\\{m<r\leq n:\Theta(r)=(x,y,z)\\}\end{split}$ (9)

Then $\bar{S}_{(0,0,1)}=\bar{S}$ and $\bar{S}_{(0,0,1)}(n)=\bar{S}(n)$.

Note that our set $\bar{S}_{(x,y,z)}$ is equivalent to the set
$H^{(\alpha,\beta)}$ used by Hajdu and Papp, where $\alpha=x$ and
$\beta=3^{y}(-1)^{z}$.

We note that the results from the previous section show that

$\bar{S}(2^{k})=1/8*2^{k}+\mathcal{O}(\sqrt{2^{k}})\text{ \, \, and \, \,
}\bar{S}(3*2^{k})=1/8(3*2^{k})+\mathcal{O}(\sqrt{2^{k}})$

The same bound can be shown to apply for other sets of numbers defined by
$(\gamma,\alpha_{3},\alpha_{5})$. In summary, as for the case when
$(x,y,z)=(0,0,1)$, the set $\bar{S}_{(x,y,z)}(n)$ is accepted by an automaton
constructed from the automata "gamma", "a3" and "a5" displayed earlier. For
example, the set of $\\{n:\Theta(n)=(1,1,1)\\}$ is accepted by the automaton
defined by the Walnut command:

    def factauto111 "?lsd_2 ~$gamma(n) & ~$a3(n) & ~$a5(n)":

As before, a linear representation of this automaton in terms of the matrices
$M_{0}$ and $M_{1}$ and the vectors $v$ and $w$ comes from the command:

    eval sumfact111 n "?lsd_2 (j>=1) & (j<=n) & $factauto111(j)":

For all choices of $(x,y,z)$, the matrix $M_{0}$ that Walnut produces is the
same as the matrix $M_{0}$ that was produced by the automaton "sumfact" that
was used in analysing the $\Theta(n)=(0,0,1)$ case. All the previous analysis
related to the case $\Theta(n)=(0,0,1)$ therefore applies to the other choices
of $(x,y,z)$. As a result, the same methods show that, for any $(x,y,z)$,

$\bar{S}_{(x,y,z)}(2^{k})=1/8*2^{k}+\mathcal{O}(\sqrt{2^{k}})\text{ \, \, and
\, \, }\bar{S}_{(x,y,z)}(3*2^{k})=1/8(3*2^{k})+\mathcal{O}(\sqrt{2^{k}}).$
(10)

It then follows that, for any choice of $(x,y,z)$ and $k$,

$\begin{split}\bar{S}_{(x,y,z)}(2^{k},2^{k+1})&=1/8*2^{k}+\mathcal{O}(\sqrt{2^{k}})\\\
\bar{S}_{(x,y,z)}(2^{k+1},3*2^{k})&=1/8*2^{k}+\mathcal{O}(\sqrt{2^{k}})\\\
\bar{S}_{(x,y,z)}(3*2^{k-2},2^{k})&=1/8*2^{k}+\mathcal{O}(\sqrt{2^{k}}).\end{split}$
(11)

######  Lemma 5.1 (Hajdu and Papp Lemma 3.1).

For any integer $k$ with $k\geq 1$, we have
$\gamma(t*2^{k}+i)=\gamma(i)+\gamma(t*2^{k})\pmod{2}$ for $0\leq i<2^{k}$.

For the following lemma, we divide the region $[0,2^{s})$ into four equal
sized subregions given by:

$\begin{split}I_{1}:=&\\{i:0\leq i<2^{s-2}\\}\\\ I_{2}:=&\\{i:2^{s-2}\leq
i<2^{s-1}\\}\\\ I_{3}:=&\\{i:2^{s-1}\leq i<3*2^{s-2}\\}\\\
I_{4}:=&\\{i:3*2^{s-2}\leq i<2^{s}\\}.\end{split}$ (12)

######  Lemma 5.2 (Hajdu and Papp Lemma 3.3).

Let $k$ be an integer and $t$ an odd integer with $t\geq 1$. Then, for each
$j\in\\{1,2,3,4\\}$, there exist numbers $c(t,j)\in\\{1,3,5,7\\}$ such that

$3^{\alpha_{3}(t*2^{k}+i)}*(-1)^{\alpha_{5}(t*2^{k}+i)}=c(t,j)*3^{\alpha_{3}(i)}*(-1)^{\alpha_{5}(i)}\pmod{8}$

for all $i\in I_{j}$.

######  Lemma 5.3.

For integers $r,s$ and $t$ satisfying $0\leq r<s$, $0\leq t<2^{s}$, $2\not|t$,

$\bar{S}(t*2^{s},\,\,t*2^{s}+2^{r})=1/8*2^{r}+\mathcal{O}(\sqrt{2^{r}})$

###### Proof.

$\begin{split}\bar{S}(t*2^{s},\,\,t*2^{s}+2^{r})&=\\#\\{i:0\leq
i<2^{r},\Theta(t*2^{s}+i)=(0,0,1)\\}\\\ &=\sum_{j=1}^{4}\\{i\in
I_{j}:\Theta(t*2^{s}+i)=(0,0,1)\\}\end{split}$

where $I_{j}$ are the regions defined in ( 12 ). From lemma 5.1, if $0\leq
i<2^{s}$,

$\Theta_{1}(t*2^{s}+i)=0\iff\Theta_{1}(i)=-\Theta_{1}(t*2^{s})$

From lemma 5.2, if $i\in I_{j}$,

$\begin{split}(\Theta_{2}(t*2^{s}+i),\,\,\,&\Theta_{3}(t*2^{s}+i))=(0,1)\\\
&\iff c(t,j)*3^{\alpha_{3}(i)}*(-1)^{\alpha_{5}(i)}\equiv-1\pmod{8}\\\
&\iff(\Theta_{2}(i),\Theta_{3}(i))=(d_{2}(t,j),d_{3}(t,j))\end{split}$

for some $d_{1},d_{2}\in\\{0,1\\}$. Hence,

$\begin{split}\bar{S}(t*2^{s},\,\,t*2^{s}+2^{r})=&\sum_{j=1}^{4}\\#\\{i\in
I_{j}:\Theta(i)=(-\Theta_{1}(t*2^{k}),d_{2}(t,j),d_{3}(t,j))\\}\\\
=&\,\,\bar{S}_{(x,y,z)_{1}}(2^{r-2})+\bar{S}_{(x,y,z)_{2}}(2^{r-2},2^{r-1})\,\,+\\\
&+\bar{S}_{(x,y,z)_{3}}(2^{r-1},3*2^{r-2})+\bar{S}_{(x,y,z)_{4}}(3*2^{r-2},2^{r})\\\
=&\,\,4*(1/8*2^{r-2}+\mathcal{O}(\sqrt{2^{r}}))\\\
=&\,\,1/8*2^{r}+\mathcal{O}(\sqrt{2^{r}}).\end{split}$

where $(x,y,z)_{j}=(-\Theta_{1}(t*2^{k}),d_{2}(t,j),d_{3}(t,j))$ for
$j\in\\{1,2,3,4\\}$ and we have used (11). ∎

The argument in the following theorem comes from theorem 2.3 in the paper by
Hajdu and Papp.

######  Theorem 5.4.

$\bar{S}(n)=1/8*n+\mathcal{O}(\sqrt{n})$

###### Proof.

We first note that (10) implies that for any $(x,y,z)$,

$\begin{split}\bar{S}_{(x,y,z)}(2^{k-2},2^{k-1})=&1/8*2^{k-2}+\mathcal{O}(\sqrt{2^{k}})\\\
\bar{S}_{(x,y,z)}(2^{k-1},3*2^{k-2})=&1/8*2^{k-2}+\mathcal{O}(\sqrt{2^{k}})\\\
\bar{S}_{(x,y,z)}(3*2^{k-2},2^{k})=&1/8*2^{k-2}+\mathcal{O}(\sqrt{2^{k}}).\end{split}$
(13)

Let $n$ be an integer with binary representation $n=\sum_{1}^{j}2^{f_{i}}$,
where $f_{1}>\dots>f_{j}\geq 0$. Then,

$\bar{S}(n)=\bar{S}(2^{f_{1}})+\bar{S}(2^{f_{1}},2^{f_{1}}+2^{f_{2}})+\dots+\bar{S}(2^{f_{1}}+\dots+2^{f_{j-1}},2^{f_{1}}+\dots+2^{f_{j}})$

Any of the terms above can be written as

$\bar{S}(2^{f_{1}}+\dots+2^{f_{l-1}},\,\,2^{f_{1}}+\dots+2^{f_{j}})=\bar{S}(t*2^{f_{l}},\,\,t*2^{f_{l}}+2^{f_{l+1}})$

where $t<2^{f_{l}}$ is odd. Since,

$\bar{S}(t*2^{f_{l}},\,\,t*2^{f_{l}}+2^{f_{l+1}})=1/8*2^{f_{l+1}}+\mathcal{O}(\sqrt{2}^{f_{l}+1})$

by lemma 5.3, we have

$\begin{split}\bar{S}(n)=&\sum_{l=1}^{j}1/8*2^{f_{l}}+\mathcal{O}(\sqrt{2}^{f_{l}})\\\
=&1/8*n+\mathcal{O}((\sqrt{2}^{f_{1}+1}-1)/(\sqrt{2}-1))\\\
=&1/8*n+\mathcal{O}(\sqrt{2}^{f_{1}}).\end{split}$

Since $f_{1}\leq\log_{2}(n)$ we have

$\bar{S}(n)=1/8*n+\mathcal{O}(\sqrt{n}).$

∎

## References

  * [1] Rob Burns.  Factorials and Legendre’s three-square theorem.  arXiv, 2021.
  * [2] Jean-Marc Deshouillers and Florian Luca.  How often is n! a sum of three squares?  In The Legacy of Alladi Ramakrishnan in the Mathematical Sciences, pages 243–251. Springer New York, 2010.
  * [3] Chen Fei Du, Hamoon Mousavi, Eric Rowland, Luke Schaeffer, and Jeffrey Shallit.  Decision algorithms for fibonacci-automatic words, ii: Related sequences and avoidability.  Theoretical Computer Science, 657:146–162, 2017.
  * [4] Chen Fei Du, Hamoon Mousavi, Luke Schaeffer, and Jeffrey Shallit.  Decision algorithms for fibonacci-automatic words, III: Enumeration and abelian properties.  International Journal of Foundations of Computer Science, 27(08):943–963, dec 2016.
  * [5] Daniel Goč, Hamoon Mousavi, and Jeffrey Shallit.  On the number of unbordered factors.  In Language and Automata Theory and Applications, pages 299–310. Springer Berlin Heidelberg, 2013.
  * [6] Andrew Granville and Yiliang Zhu.  Representing binomial coefficients as sums of squares.  The American Mathematical Monthly, 97(6):486–493, 1990.
  * [7] L. Hajdu and A. Papp.  On asymptotic density properties of the sequence n!  Acta Arithmetica, 184(4):317–340, 2018.
  * [8] Hamoon Mousavi.  Automatic theorem proving in Walnut.  arXiv, 2016.
  * [9] Hamoon Mousavi, Luke Schaeffer, and Jeffrey Shallit.  Decision algorithms for fibonacci-automatic words, i: Basic results.  RAIRO Theor. Informatics Appl., 50:39–66, 2016.
  * [10] Narad Rampersad and Jeffrey Shallit.  Congruence properties of combinatorial sequences via walnut and the rowland-yassawi-zeilberger automaton.  arXiv, 2021.
  * [11] Jeffrey Shallit.  Sumsets of wythoff sequences, fibonacci representation, and beyond.  arXiv, 2020.
  * [12] Jeffrey Shallit.  Frobenius numbers and automatic sequences.  arXiv, 2021.
  * [13] Jeffrey Shallit.  Intertwining of complementary thue-morse factors.  arXiv, 2022.
  * [14] The Sage Developers.  SageMath, the Sage Mathematics Software System (Version 8.2), 2018.  https://www.sagemath.org.