# Synthetic magnetism for solitons in optomechanical array

P. Djorwé<EMAIL_ADDRESS>Department of Physics, Faculty of Science,
University of Ngaoundere, P.O. Box 454, Ngaoundere, Cameroon  Stellenbosch
Institute for Advanced Study (STIAS), Wallenberg Research Centre at
Stellenbosch University, Stellenbosch 7600, South Africa  H. Alphonse
<EMAIL_ADDRESS>Department of Marine Engineering, Limbé Nautical Arts and
Fisheries Institute, P.O Box 485, Limbé, Cameroon  S. Abbagari
<EMAIL_ADDRESS>Department of Basic Science, National Advanced School of
Mines and Petroleum Industries, The University of Maroua, P.O Box 08, kaélé,
Cameroon  S. Y. Doka<EMAIL_ADDRESS>Department of Physics, Faculty of
Science, University of Ngaoundere, P.O. Box 454, Ngaoundere, Cameroon  S. G.
Nana Engo<EMAIL_ADDRESS>Department of Physics, Faculty
of Science, University of Yaounde I, P.O. Box 812, Yaounde, Cameroon

###### Abstract

We propose a synthetic magnetism to generate and to control solitonic waves in
$\rm{1D}-$optomechanical array. Each optomechanical cavity in the array
couples to its neighbors through photon and phonon coupling. We create the
synthetic magnetism by modulating the phonon hopping rate through a modulation
frequency, and a modulation phase between resonators at different sites. When
the synthetic magnetism effect is not taken into account, the mechanical
coupling play a scrucial role of controlling and switching the waves from
bright to dark solitons, and it even induces rogue wave-like a shape in the
array. For enough mechanical coupling strength, the system enters in a strong
coupling regime through splitting/crossing of solitonic waves leading to
multiple waves propagation in the array. Under the synthetic magnetism effect,
the phase of the modulation enables a good control of the wave propagation,
and it also switches soliton shape from bright to dark, and even induces rogue
waves as well. Similarly to the mechanical coupling, the synthetic magnetism
offers another flexible way to generate plethora of solitonic waves for
specific purposes. This work opens new avenues for optomechanical platforms
and sheets light on their potentiality of controlling and switching solitonic
waves based on synthetic magnetism.

Optomechanics, Solitons, synthetic magnetism

###### pacs:

42.50.Wk, 42.50.Lc, 05.45.Xt, 05.45.Gg

##  I Introduction

Optomechanical cavities (COM) have been revealed as being a nice platform for
solitonic waves investigation owing to their intrinsic nonlinear interactions
Leijssen et al. (2017); Doolin et al. (2014); Djorwe et al. (2019a); Kingni et
al. (2020); Djorwé et al. (2013a). Nonlinear interactions in simple COM have
led to numerous interesting studies such as quantum entanglement Djorwé et al.
(2014); Kotler et al. (2021), squeezing Djorwé et al. (2013b); Wollman et al.
(2015), high-order sideband combs generation Cao et al. (2016); Djorwe et al.
(2020); Miri et al. (2018), metamorphoses of basin boundaries Djorwe et al.
(2022), and chaos Monifi et al. (2016). To harness more features related to
optomechanical coupling, simple COMs have been extended to optomechanical
array, and this has led to exotic physics including synchronization of COMs
Djorwé et al. (2020); Sheng et al. (2020), Parity-Time symmetry El-Ganainy et
al. (2018), sensors based on exceptional points Hodaei et al. (2017); Chen et
al. (2017); Djorwe et al. (2019b), solitonic waves Alphonse et al. (2022), and
bosonic synthetic magnetism for topological transport phenomena Fang et al.
(2012); Schmidt et al. (2015); Fang et al. (2017); Brendel et al. (2017); Wang
et al. (2020).

Unlike electrons, photons and phonons do not feel a magnetic field, because
they are not charged. As such, most of the physics connected with charged
particles in magnetic fields is missing for bosons. To feed this gap,
researchers started to engineer artificial magnetic fields for bosons to make
them akin to electrons in a magnetic field. Different schemes for synthetic
magnetism have been proposed in the literature for specific purposes. These
proposals are mostly related to array of COMs, where cells are connected
through photon and/or phonon coupling. The engineering of the synthetic field
lies on the modulation of this photon/phonon hopping rate. In Schmidt et al.
(2015), an artificial magnetic fields for photons has been created in
optomechanical crystal. This magnetic field results from a modulation of
photon coupling on the lattice, and was used to achieve photon transport in
the presence of an artificial Lorentz force, edge states, and the photonic
Aharonov–Bohm effect. Pseudomagnetic field for sound has been created in
Brendel et al. (2017) for chiral transport of sound waves at nanoscale. The
resulted Pseudomagnetic field was created by simple geometrical modifications
of a snowflake phononic crystal. Along the same line, a magnetic gauge field
for phonons has been demonstrated in a sliced optomechanical crystal nanobeam
in Mathew et al. (2020). The gauge field were induced through multi-mode
optomechanical interactions, and has been used to break time-reversal
symmetry. Under a dynamical modulation of the system, the authors observed
nonreciprocal phonon transport between resonators of different frequencies
mediated by radiation pressure forces. In Wang et al. (2020), a novel
mechanism to generate synthetic magnetic field for phonon lattice by Floquet
engineering auxiliary qubits has been proposed to realize circulator transport
and analog Aharonov–Bohm effects for acoustic waves. Owing to these
interesting physics fostered by synthetic magnetism in optomechanics, here we
propose to investigation solitonic waves in 1D optomechanical array which
hosts both photon and phonon coupling between cells. In our proposal, the
phonon coupling between cells is harmonically modulated in time, and this
induces artificial magnetic field for phonons. Such a mechanism has been
realized in Fang et al. (2012) to create a effective magnetic field for
photons by controlling the phase of the dynamical modulation. By tunning the
phase of the modulation, we observed i) a perfect control and switching
between bright and dark mechanical solitons along the array, ii) a strong
coupling regime inducing rogue waves, and iii) multiple solitons
crossings/splittings. These interesting features reveal the scrucial effect of
the synthetic magnetism on the solitonic waves, and can be extended to similar
systems. This work opens new avenues for solitons generation, and provides a
smoothly way to switch from bright to dark solitons depending on the purposes.

The rest of this work is organized as follows. In section II, we describe the
system and the related equations. The mechanical coupling and the synthetic
magnetism effect on the solitonic waves are presented in section III and we
conclude the work in section IV.

##  II Model and dynamical equations

Our proposal consists on $N$ cells of COM which are both optically and
mechanically coupled to each other as sketched in Figure 1.

Figure 1: Dynamically modulated $\rm{1D}-$optomechanical array for an
effective magnetic field for phonons. Each cell supports a mechanical
resonator (blue color) and an optical cavity (red color). The neighbors’ cell
are optically (mechanically) coupled through the hopping coupling rate $J_{0}$
($J_{m}$). The mechanically hopping rate is dynamically modulated as
$J_{m}(\phi)=J_{m}\cos{(\Omega t+\phi)}$.

The mechanical resonator and the optical cavity of each cell are depicted by
the blue and red color respectively, and the photon (phonon) hopping coupling
rate between neighboring cells is captured by $J_{0}$ ($J_{m}$). The dynamics
of such an array of COM is captured by the Hamiltonian (with $\hslash=1$),

$H=H_{0}+H_{J_{0}}+H_{J_{m}}+H_{int},$ (1)

where

$\displaystyle H_{0}$
$\displaystyle=\sum_{n}\left(\omega_{m,n}b^{\dagger}_{n}b_{n}+\omega_{c,n}a^{\dagger}_{n}a_{n}\right),$
(2) $\displaystyle H_{J_{0}}$
$\displaystyle=-\frac{J_{0}}{2}\sum_{n}\left(a^{\dagger}_{n}a_{n+1}+a_{n}a^{\dagger}_{n+1}\right),$
(3) $\displaystyle H_{J_{m}}$
$\displaystyle=-\frac{J_{m}(\phi)}{2}\sum_{n}\left(b^{\dagger}_{n}b_{n+1}+b_{n}b^{\dagger}_{n+1}\right),$
(4) $\displaystyle H_{int}$
$\displaystyle=-\sum_{n}g_{n}a^{\dagger}_{n}a_{n}\left(b^{\dagger}_{n}+b_{n}\right),$
(5)

where $H_{0}$, $H_{J_{0}}(H_{J_{m}})$, and $H_{int}$ are the free
optomechanical, the photon (phonon) hopping and optomechanical interaction
Hamiltonian, respectively. The optical and mechanical mode are described by
their $n^{th}$ related bosonic field operators, $a_{n}$ ($a^{\dagger}_{n}$)
and $b_{n}$ ($b^{\dagger}_{n}$), respectively. The mechanical displacement is
connected to its operators through $x_{n}=x^{n}_{ZPF}(b_{n}+b^{\dagger}_{n})$
where $x^{n}_{ZPF}$ is the zero-point fluctuation amplitude of the $n^{th}$
mechanical resonator. The other parameters are the mechanical (cavity)
frequency $\omega_{m,n}$ ($\omega_{c,n}$), and the optomechanical coupling
$g_{n}$. The term $J_{m}(\phi)$ in Equation 4 is the modulated phonon hopping
term, which is introduced to reveal the synthetic magnetism features. This
term is explicitly given by $J_{m}(\phi)=J_{m}\cos{(\Omega t+\phi)}$, where
$\Omega$ is the modulation frequency, and $\phi$ is the phase of the
modulation between resonators at different sites. From now on, we assume that
our cells are degenerated, and therefore use a simplified notation by dropping
the index $n$ such as $\omega_{m,n}\equiv\omega_{m}$,
$\omega_{c,n}\equiv\omega_{c}$, and $g_{n}\equiv g$. From the Hamiltonian
given in Equation 1, one deduces the following Quantum Langevin Equations
(QLEs),

$\begin{split}\dot{a}_{n}&=i\left[-\omega_{c}+g(b^{\dagger}_{n}+b_{n})\right]a_{n}+i\frac{J_{0}}{2}(a_{n+1}+a_{n-1}),\\\
\dot{b}_{n}&=-i\omega_{m}b_{n}+ig|a_{n}|^{2}+i\frac{J_{m}(\phi)}{2}(b_{n+1}+b_{n-1}).\end{split}$
(6)

To capture the soliton dynamics from Equation 6, we assume that photon and
phonon numbers are large enough to consider classical limit where noise terms
can be neglected. Moreover, the operators can be substituted by their $c$-mean
values $\langle a_{n}\rangle=\alpha_{n}$ and $\langle b_{n}\rangle=\beta_{n}$.
In the rotating frame moving at $\omega_{c}+J_{0}$, the above QLEs yields,

$\begin{split}\dot{\alpha}_{n}&=i[g(\beta^{\dagger}_{n}+\beta_{n})]\alpha_{n}+i\frac{J_{0}}{2}(\alpha_{n+1}+\alpha_{n-1}-2\alpha_{n}),\\\
\dot{\beta}_{n}&=-i\omega_{m}\beta_{n}+ig|\alpha_{n}|^{2}+i\frac{J_{m}(\phi)}{2}(\beta_{n+1}+\beta_{n-1}-2\beta_{n}),\end{split}$
(7)

which is the set of equations that will be considered latter on in our
numerical simulations.

##  III Soliton dynamics under the mechanical coupling and the synthetic
magnetism

To carry out numerical analysis of our proposal, we consider the dynamical
equations given in Equation 7. We seek to harness the synthetic magnetic field
effect on the solitonic wave propogation under the initial solution,

$\alpha_{n}(0)=\tfrac{\sqrt{\frac{\omega_{m}J_{0}}{2}}}{gx_{0}}\operatorname{sech}\left(\frac{n-x(0)}{x_{0}}\right)\exp(i\rm{k}n),$
(8)

where $x_{0}$ and $x(0)$ are the soliton width and the initial position of the
envelope respectively, and $\rm{k}$ is the wave vector. Throughout the
numerical simulations, the initial conditions $\alpha_{n}=\alpha_{n}(0)$ and
$\beta_{n}(0)=0$ will be used for the optical field and the mechanical
resonator, respectively. Since we are interested on the mechanical coupling
and synthetic magnetism effects, the wave vector $\rm{k}$ will be set to zero
( $\rm{k=0}$).

###  III.1 Soliton dynamics under the mechanical coupling effects

We are interesting here in analyzing the effects of the mechanical coupling on
solitonic wawes. Therefore, the synthetic magnetism effect is not acconted
($\Omega=0;\phi=0$). Figure 2 shows the $\rm{2D}$ projection profile of
mechanical soliton emerging in our system. These figures correspond to
different values of the mechanical coupling $J_{m}$ when the effect of the
synthetic magnetism is not taken into account. The values of $J_{m}$ for this
set of figures are $\tfrac{J_{m}}{\omega_{m}}=0;0.5;0.7$ and $0.8$ for Figure
2 (a;b;c and d), respectively. Different shapes of the mechanical soliton can
be seen, and all are propagating at the center of the array. It can be
observed that there are more and stable solitonic harmonics propagating in
Figure 2a than in Figure 2b.

Figure 2: Dynamics of bright solitons in index and time $\rm{2D}$
representation. The mechanical coupling from (a) to (d) are
$\tfrac{J_{m}}{\omega_{m}}=0;0.5;0.7$ and $0.8$, respectively. The rest of the
parameters are $N=100$, $J_{0}=2\omega_{m}$, $g=10^{-4}\omega_{m}$,
$x_{0}=10$, $n=100$, $x(0)=\frac{n}{2}$, $\Omega=0$, $\phi=0$ and $k=0$.

In Figure 2c and Figure 2d, the soliton peaks are less present on the array,
and their high is also decreasing compared to those in Figure 2a and Figure
2b. More importantly, one observes that some hollows/deeps are appearing on
Figure 2c and Figure 2d. These hollows reveal the fact that there are new
solitons which are propagating with reversed peaks in the array, and such a
feature is reminiscent of dark solitons. This transition form bright to dark
solitonic waves can be further pointed out by keeping increase the mechanical
coupling as depicted in Figure 3. From Figure 3a to Figure 3d, the phonon
hopping rate is ramping from $\tfrac{J_{m}}{\omega_{m}}=1.01;1.1;1.3$ to
$1.4$.

Figure 3: Dynamics of dark solitons in index and time $\rm{2D}$
representation. The mechanical coupling from (a) to (d) are
$\tfrac{J_{m}}{\omega_{m}}=1.01;1.1;1.3$ and $1.4$, respectively. The other
parameters are the same as in Figure 2.

By comparing Figure 2 and Figure 3, one observe that the background of these
set of figures has switched from dark to bright, meaning that the peaks on
these two configurations are globally reversed one to another, and this
reveals that the transitional phase from bright to dark solitonic waves has
been realized. In Figure 3a, it can be seen that the increase of $J_{m}$ has
led to the generation of two separated dark solitons. As this coupling keep
increase, it induces the crossing/splitting of these two dark solitons as
shown in Figure 3b, which means that the system has reached a strong coupling
regime. For strong enough strength of $J_{m}$, the array exhibits complex
features as depicted in Figure 3c and Figure 3d. These complex features seem
to exhibit a combination of both bright and dark solitons, which is a
signature of rogue waves. To highlight the transitional phase between
bright/dark solitonic wave and the rogue wave features, we have exemplified
some aspects of Figure 3a and Figure 3b in Figure 4. Figure 4a is a
respresentation of Figure 3a at the longest time considered here
($t=t_{max}=50\omega_{m}$). This figure shows how the maximum peack of the
soliton is around zero, which agrees with the observation made in Figure 2
pointing out that an increase of the mechanical coupling $J_{m}$ induces a
decrease of the soliton’s peack. Moreover, the two hollows/deeps which started
from Figure 2c (and amplified in Figure 3a) are clearly illustrated in Figure
4a which shows how the solitonic wave is at the transitional threshold between
bright to dark solitons. Indeed, increasing $J_{m}$ will totally lead to
negative soliton’s peack. We conclude that this phase transition happens at
the vicinity of $J_{m}=\omega_{m}$. Figure 4b is a $\rm{3D}$-representation of
Figure 3d (for $J_{m}=1.4\omega_{m}$), which is well above the transitional
threshold. As predicted, this figure depicts a mixing of bright and dark
solitons which is a rogue waves signature. Owing to these features, we
conclude that the mechanical coupling strength can be used to control
solitonic waves propagation and to switch from bright to dark solitons in
optomechanical array.

Figure 4: (a) is the representation of Figure 3a at $t=50\omega_{m}$ to
highlight the phase transition between bright and dark soliton. (b) is a
$\rm{3D}$ illustration of Figure 3d to depict rogue waves feature in strong
coupling regime ($J_{m}=1.4\omega_{m}$).

To point out the optomechanical nature of our study, we would like to carry
out the effect of the mechanical coupling on the optical solitons as well.
This analysis is shown in Figure 5 where the optical coupling strength is
maintained at $J_{0}=2\omega_{m}$. As it can be seen, there is a stable
optical solitonic wave that is propagating at the center of the array in
Figure 5a, which corresponds to the case without mechanical coupling
$J_{m}=0$. The effect of the phonon hopping can be see from Figure 5b, where
the width of the wave is shrinking and it becomes blurred in Figure 5c as the
phonon hopping rate increases from $\tfrac{J_{m}}{\omega_{m}}=0.8$ to $1.01$.
This blurred shape of the optical solitons is the sign that their peak is
decreasing, meaning that the wave is lossing its energy that is spread to
cells located far from the array’s center. To confirm such statement, we have
further increases phonon hopping strength ($\tfrac{J_{m}}{\omega_{m}}=1.1$),
and we observed an optical soliton splitting in Figure 5c. This splitting
reveals that our system has reached a strong coupling regime, where solitons
can cross or split leading to multiple solitonic waves propagating in the
array. Indeed, Figure 5(c) to Figure 5e are plotted respectively for
$\tfrac{J_{m}}{\omega_{m}}=1.1;1.2;1.3$, which clearly are above the strong
coupling regime, and we can see two (see Figure 5c), three (see Figure 5d) and
four (see Figure 5e) solitonic waves propagating in the array. From this
analysis, we conclude that the strong coupling induced by the mechanical
coupling is reflected both in our optical and mechanical resonators through
the splitting/crossing of the solitonc waves. This interplay between optical
and mechanical waves in our array is mediated by the optomechanical coupling
$g$, which is the key parameter that connects optical and mechanical
resonators in each cells of the array.

Figure 5: Dynamics of the optical solitonic waves showing solitonic waves
splittings/crossings due to a strong coupling. From (a) to (f), the mechanical
coupling is respectively, $\tfrac{J_{m}}{\omega_{m}}=0;0.8;1.01$, $1.1$, $1.2$
and $1.3$. The other parameters are the same as in Figure 2. One can see two
(c), three (d) and four (e) solitonic waves crossing/splitting in the array.

###  III.2 Soliton dynamics under the synthetic magnetism effect

To point out the effect of the synthetic magnetism, we can either set on the
modulation frequency $\Omega$ or the phase of the modulation between
resonators at different sites $\phi$. For a proof of concept purpose and
without lost of generality, we will not set on these two parameters
simultaneously in this work since they have similar effects. However, we will
mostly focus on the phase-induced effects. Moreover, we would like to point
out that i) for $(\Omega t+\phi)\equiv\frac{n\pi}{2}$ with n
$\in\mathbb{Z^{\ast}}$, our system only hosts the optical coupling between
cells and the related results have been shown in Alphonse et al. (2022); ii)
for $(\Omega t+\phi)\equiv n\pi$ with n $\in\mathbb{Z}$, our proposal does not
host synthetic magnetism and the related results reduce to those shown in
Figure 2 and Figure 3. The effect of the synthetic magnetic field is depicted
in Figure 6, where the mechanical coupling is fixed to
$\tfrac{J_{m}}{\omega_{m}}=0.5$ for Figure 6(a-c) and to
$\tfrac{J_{m}}{\omega_{m}}=1.5$ for Figure 6(d-e). To point out the effect of
the modulation frequency, we have tunned it to $\Omega=\tfrac{2\pi\nu}{n}$
where the frequency is $\nu=50\rm{Hz}$, with $n\in\mathbb{Z^{\ast}}$. Starting
from Figure 2b and by setting $n=3$ and $\phi=0$, we obtained Figure 6a. It
can be seen in Figure 6a that the peaks have spread along the center of the
array, meaning that $\Omega$ enhances propagation of solitonic waves. Similar
effect is observed in Figure 6b, where we have set $\Omega=0$ and $\phi=\pi$.
As we previously predicted, Figure 6c depicts the solitonic wave propagation
when the effect of the mechanical coupling is annihilated for
$\phi=\frac{\pi}{2}$, leading to a well known result (see Figure 2a). To
further figure out the effect of the phase ($\phi$) on our system, let us
still start with a stable configuration corresponding to $\Omega=0$, $\phi=0$
and $J_{m}=1.5\omega_{m}$ (not shown). From such a stable solitonic wave, we
tune the phase to $\phi=\frac{\pi}{3};\frac{\pi}{4}$ and $\frac{\pi}{6}$. The
corresponding outputs are respectively depicted on Figure 6(d-e), and it be
clearly seen that like $J_{m}$, the phase induces a transition from bright to
dark solitons. Furthermore, the rogue wave-like feature can be also recovered
through phase adjustment. We can thus conclude that similar to the phonon
hopping effect, the phase $\phi$ is equaly an important parameter to control
solitonic waves and to smoothly switch from bright to dark soliton features.

Figure 6: Synthetic magnetism effect on the soliton dynamics. (a) corresponds
to $\Omega=\tfrac{2\pi\nu}{3}$ and $\phi=0$ while in (b) to (f), $\Omega=0$
and the phase is respectively fixed to
$\phi=\pi;\frac{\pi}{2};\frac{\pi}{3};\frac{\pi}{4}$ and $\frac{\pi}{6}$. The
mechanical coupling is fixed to $\tfrac{J_{m}}{\omega_{m}}=0.5$ for (a)-(c)
and to $\tfrac{J_{m}}{\omega_{m}}=1.5$ for (d)-(e). The other parameters are
the same as in Figure 2.

##  IV Conclusion

This work investigated the effect of mechanical coupling and a synthetic
magnetic field on solitonic waves propagating in $\rm{1D-}$ optomechanical
array. Within the array, each optomechanical cavity is coupled to its
neighbors through photon and phonon couplings. The synthetic magnetism is
created by modulating the phonon hopping rate. We have shown that the tunning
of the mechanical coupling leads to several solitonic waves features ranging
from bright to dark solitons, and even to the emergence of rogue waves in the
array. Moreover, this phonon hopping strength induces a strong coupling regime
revealed through splitting/crossing of solitonic waves leading to multiple
waves propagating in the array. For a fixed strength of the mechanical
coupling, we have also shown that the frequency and the phase of the synthetic
magnetism are scrucial parameters to control the solitonc wave as well. By
focusing on the phase, we were able to not only control the wave propagation,
but this phase has been proven to switch soliton shape from bright to dark,
and even to induce rogue waves in the array. This work opens new avenues for
optomechanical platforms and provides ingrediens based on phononic coupling
and synthetic magnetism to control and smoothly switch solitonic waves from
one shape to another.

## References

  * Leijssen et al. (2017) R. Leijssen, G. R. L. Gala, L. Freisem, J. T. Muhonen, and E. Verhagen, Nature Communications 8 (2017).
  * Doolin et al. (2014) C. Doolin, B. D. Hauer, P. H. Kim, A. J. R. MacDonald, H. Ramp, and J. P. Davis, Physical Review A 89 (2014).
  * Djorwe et al. (2019a) P. Djorwe, Y. Pennec, and B. Djafari-Rouhani, Scientific Reports 9 (2019a).
  * Kingni et al. (2020) S. T. Kingni, C. Tchodimou, D. P. Foulla, P. Djorwe, and S. G. N. Engo, The European Physical Journal Special Topics 229, 1117 (2020).
  * Djorwé et al. (2013a) P. Djorwé, J. H. T. Mbé, S. G. N. Engo, and P. Woafo, The European Physical Journal D 67 (2013a).
  * Djorwé et al. (2014) P. Djorwé, S. G. N. Engo, and P. Woafo, Physical Review A 90 (2014).
  * Kotler et al. (2021) S. Kotler, G. A. Peterson, E. Shojaee, F. Lecocq, K. Cicak, A. Kwiatkowski, S. Geller, S. Glancy, E. Knill, R. W. Simmonds, et al., Science 372, 622 (2021).
  * Djorwé et al. (2013b) P. Djorwé, S. N. Engo, J. T. Mbé, and P. Woafo, Physica B: Condensed Matter 422, 72 (2013b).
  * Wollman et al. (2015) E. E. Wollman, C. U. Lei, A. J. Weinstein, J. Suh, A. Kronwald, F. Marquardt, A. A. Clerk, and K. C. Schwab, Science 349, 952 (2015).
  * Cao et al. (2016) C. Cao, S.-C. Mi, T.-J. Wang, R. Zhang, and C. Wang, IEEE Journal of Quantum Electronics 52, 1 (2016).
  * Djorwe et al. (2020) P. Djorwe, J. Y. Effa, and S. G. N. Engo, Journal of the Optical Society of America B 37, A36 (2020).
  * Miri et al. (2018) M.-A. Miri, G. D’Aguanno, and A. Alù, New Journal of Physics 20, 043013 (2018).
  * Djorwe et al. (2022) P. Djorwe, J. Y. Effa, and S. G. N. Engo, Nonlinear Dynamics (2022).
  * Monifi et al. (2016) F. Monifi, J. Zhang, Ş. K. Özdemir, B. Peng, Y. xi Liu, F. Bo, F. Nori, and L. Yang, Nature Photonics 10, 399 (2016).
  * Djorwé et al. (2020) P. Djorwé, Y. Pennec, and B. Djafari-Rouhani, Physical Review B 102, 155410 (2020).
  * Sheng et al. (2020) J. Sheng, X. Wei, C. Yang, and H. Wu, Physical Review Letters 124 (2020).
  * El-Ganainy et al. (2018) R. El-Ganainy, K. G. Makris, M. Khajavikhan, Z. H. Musslimani, S. Rotter, and D. N. Christodoulides, Nature Physics 14, 11 (2018).
  * Hodaei et al. (2017) H. Hodaei, A. U. Hassan, S. Wittek, H. Garcia-Gracia, R. El-Ganainy, D. N. Christodoulides, and M. Khajavikhan, Nature 548, 187 (2017).
  * Chen et al. (2017) W. Chen, Ş. K. Özdemir, G. Zhao, J. Wiersig, and L. Yang, Nature 548, 192 (2017).
  * Djorwe et al. (2019b) P. Djorwe, Y. Pennec, and B. Djafari-Rouhani, Physical Review Applied 12, 024002 (2019b).
  * Alphonse et al. (2022) H. Alphonse, P. Djorwe, S. Abbagari, S. Y. Doka, and S. N. Engo, Chaos, Solitons and Fractals 154, 111593 (2022).
  * Fang et al. (2012) K. Fang, Z. Yu, and S. Fan, Nature Photonics 6, 782 (2012).
  * Schmidt et al. (2015) M. Schmidt, S. Kessler, V. Peano, O. Painter, and F. Marquardt, Optica 2, 635 (2015).
  * Fang et al. (2017) K. Fang, J. Luo, A. Metelmann, M. H. Matheny, F. Marquardt, A. A. Clerk, and O. Painter, Nature Physics 13, 465 (2017).
  * Brendel et al. (2017) C. Brendel, V. Peano, O. J. Painter, and F. Marquardt, Proceedings of the National Academy of Sciences 114 (2017).
  * Wang et al. (2020) X. Wang, H.-R. Li, and F.-L. Li, New Journal of Physics 22, 033037 (2020).
  * Mathew et al. (2020) J. P. Mathew, J. del Pino, and E. Verhagen, Nature Nanotechnology 15, 198 (2020).


# The complete control of scattering waves in multi-channel structures

Qi Gao  Yun-Song Zhou<EMAIL_ADDRESS>Li-Ming Zhao  Department of
Physics, Capital Normal University, Beijing, 100048. China

###### Abstract

The issue of photon spin Hall effect was generalized as a universal question
of how to control all the scattering waves in a multi-channel structure
(complete control). A general theory was proposed, which provides a simple way
to achieve the complete control. This theory shows also that the necessary
condition for complete control is that the structure must contain a complete
set of sources. To demonstrate the application of the theory, the typical
scattering patterns in the two-channel and four-channel structures are
achieved theoretically. Previous this research, one could only artificially
control the scattering waves in two channels out of a four-channel structure.

The new kind of photonic spin Hall effect (PHE) attracted great interest[1, 2,
3, 4, 5, 6, 7, 8, 9, 10] because of its application prospects and interesting
theoretical problems[11, 12, 13, 14, 15, 16, 17]. Considering this problem
from a general perspective, it belongs to the scattering control problem.
Taking the PHE on metal surface[1, 2] as an example, it is controlled by the
polarization state of incident light. The surface polariton plasmon (SPP) wave
can be controlled traveling to left (left channel) or right (right channel)
artificially, here a channel corresponds to a guide mode. So, the metal
surface is a two-channel structure, the two channels are separated at the
point where the scatterer is located. Most of the works currently are focused
on the two-channel structures, but some researches are also extended to the
four-channel structures[16, 17, 18, 19, 20]. However, if we review these works
on four-channel structures, a common point can be found that, only two
channels can be selected from the four channels, in which scattering waves are
controllable as people whish. The other waves in the remaining two channels
can not be controlled artificially (known later that the waves in remaining
channels are determined by the waves in the selected channels). Of course, the
two controllable channels are optional. This kind of control can be name as
2/4 control. (noticing that the PHE on metal surface is 2/2 control because
the waves in the two channels can be distributed with any proportion). Thus, a
series of interesting questions arose: how to realize 3/4, and 4/4 controls in
a four-channel structure? The more universal questions are proposed: in a
n-channel structure (if a channel involves two modes it should be regard as
two channels), how to realize the complete control(the waves in all channels
are controllable)? and what’s the condition for achieving the complete
control?

To answer these questions, let’s first consider the simplest structure[1, 2].
A notch cut on a metal surface which will divide the surface into two
channels, channel 1 on the left side and channel 2 on the right side of the
notch. When the notch is illuminated by an incident light, two SPP waves will
be excited and transmit respectively in both channels[2]. The magnetic fields
of the eigenmodes in channel 1 and channel 2 are expressed with
$\mathbf{h}_{1}$ and $\mathbf{h}_{2}$ respectively(see Fig. 1(a)). When the
incident light is circularly polarized, the illuminated notch can be thought
as a rotational dipole[1, 2, 16, 17], which is composed of two independent
dipoles along two perpendicular directions, just as shown in Fig. 1(a). Each
independent dipole can be called as a light source. The contributions from the
two sources can be discussed independently.

For the $\mathbf{d}_{x}$ dipole (source 1), its radiation magnetic field is
$\mathbf{H}^{(\text{i})}_{1}$, which will excite
$\mathbf{H}^{(\text{s})}_{1}$, the magnetic field of SPP wave (Fig. 1(b)).
$\mathbf{H}^{(\text{s})}_{1}$ is the total scattering wave involving the both
waves in the different directions. In other words
$\mathbf{H}^{(\text{s})}_{1}=s_{11}\mathbf{h}_{1}+s_{21}\mathbf{h}_{2}$. As
the same reason, the $\mathbf{d}_{z}$ dipole (source 2) has the relationship
$\mathbf{H}^{(\text{s})}_{2}=s_{12}\mathbf{h}_{1}+s_{22}\mathbf{h}_{2}$ (Fig.
1(c)). Thus, we have the matrix equation

$\begin{aligned}
\begin{pmatrix}{{\mathbf{H}}_{1}^{({\text{s}})}}&{{\mathbf{H}}_{2}^{({\text{s}})}}\end{pmatrix}=\begin{pmatrix}{{\mathbf{h}_{1}}}&{{\mathbf{h}_{2}}}\end{pmatrix}\begin{pmatrix}{{s_{11}}}&{{s_{12}}}\\\
{{s_{21}}}&{{s_{22}}}\end{pmatrix}\equiv\begin{pmatrix}{{\mathbf{h}_{1}}}&{{\mathbf{h}_{2}}}\end{pmatrix}S\end{aligned}.$
(1)

Here $S$ is named as transform matrix. If each source strength
$\mathbf{H}^{(\text{i})}_{n}$ is normalized, or standard with a fixed value,
each element of the transform matrix can be determined by theoretical
calculation, numerical simulation, or practice measurement, and then the total
field excited by the two dipoles should be composed as
$\mathbf{H}^{(\text{s})}=\mathbf{H}^{(\text{s})}_{1}+\mathbf{H}^{(\text{s})}_{2}$
(sketched in Fig. 1(d)).

However, for a certain structure, the matrix $S$ should be determined first by
the normalized or standard $\mathbf{H}^{(\text{i})}_{n}$. In actual problems,
the source strength (depend on the incident light) may not be normalized or
standard, but $A_{1}\mathbf{H}^{(\text{i})}_{1}$ and
$A_{1}\mathbf{H}^{(\text{i})}_{2}$, the total scattering field should be
$\mathbf{H}^{(\text{s})}=A_{1}\mathbf{H}^{(\text{s})}_{1}+A_{2}\mathbf{H}^{(\text{s})}_{2}$.
Combined with Eq. 1, we get

${{\mathbf{H}}^{({\text{s}})}}=\begin{pmatrix}{\mathbf{h}_{1}}&{\mathbf{h}_{2}}\end{pmatrix}\begin{pmatrix}{{C_{1}}}\\\
{{C_{2}}}\end{pmatrix},$ (2)

where

$\begin{pmatrix}{{C_{1}}}\\\
{{C_{2}}}\end{pmatrix}=S\begin{pmatrix}{{A_{1}}}\\\ {{A_{2}}}\end{pmatrix}.$
(3)

Equation 3 is a linear transform. If $S$ is an invertible matrix, then

$\begin{pmatrix}{{A_{1}}}\\\
{{A_{2}}}\end{pmatrix}={S^{-1}}\begin{pmatrix}{{C_{1}}}\\\
{{C_{2}}}\end{pmatrix},$ (4)

Eqs. 2 and 4 provide the way to artificially control the scattering field. For
instance, if we hope the scattered light transmitting only in $x$ direction
(Fig. 1(d)), the wavefunction should be
$\mathbf{H}^{(\text{s})}=C_{1}\mathbf{h}_{1}+C_{2}\mathbf{h}_{2}$, it means
$C_{1}=0$ and $C_{2}=1$. The corresponding source strengths $A_{1},A_{2}$ can
be obtained from Eq. (4). In general case, any wanted scattering field, which
is represented by $C_{1},\ C_{2}$ can be obtained when we set the source
strength $A_{1},\ A_{2}$ according to Eq. 4. Thus, the complete control is
achieved. In other words, the waves in two channels can be controlled in the
two-channel structure, it is named as 2/2 control.

Figure 1: The photon spin Hall effect on metal surface can be thought as the
scattering control process in two-channel structure. (a) the dashed circle
denotes the rotational dipole, which can be separated with both the level and
vertical dipoles. The surface on the left (right) side of the dipole is
channel 1 (channel 2), in which the SPP mode is $\mathbf{h}_{1}$
($\mathbf{h}_{2}$); (b) the SPP wave $\mathbf{H}_{1}^{({\text{s}})}$ excited
by the level dipole. (c) the SPP wave $\mathbf{H}_{2}^{({\text{s}})}$ excited
by the vertical dipole. (d) the total wave
$\mathbf{H}^{({\text{s}})}=\mathbf{H}_{1}^{({\text{s}})}+\mathbf{H}_{2}^{({\text{s}})}$
scattered by the rotational dipole.

Now, it’s easy to extend this theory to the multi-channel structure with the
same way. For a $n$-channel structure, its $n$ eigenfunctions are
$\mathbf{h}_{1},\dots,\mathbf{h}_{n}$. If we set $n$ independent light sources
in the structure, the normalized or standard fields are expressed as
$\mathbf{H}^{(\text{i})}_{1},\dots,\mathbf{H}^{(\text{i})}_{n}$, which will
excite the scattering fields
$\mathbf{H}^{(\text{s})}_{1},\dots,\mathbf{H}^{(\text{s})}_{n}$. When the real
radiation fields are
$A_{1}\mathbf{H}^{(\text{i})}_{1},\dots,A_{n}\mathbf{H}^{(\text{i})}_{n}$, the
total scattering field will be
${{\mathbf{H}}^{({\text{s}})}}=\sum\limits_{\beta=1}^{n}{{\mathbf{H}}_{\beta}^{({\text{s}})}{A_{\beta}}}=\sum\limits_{\alpha,\beta=1}^{n}{{{\mathbf{h}}_{\alpha}}{S_{\alpha\beta}}{A_{\beta}}}=\sum\limits_{\alpha=1}^{n}{{{\mathbf{h}}_{\alpha}}{C_{\alpha}}}$.
Here

$\begin{pmatrix}{{C_{1}}}\\\ \vdots\\\
{{C_{n}}}\end{pmatrix}=\begin{pmatrix}{{s_{11}}}&\cdots&{{s_{1n}}}\\\
\vdots&\ddots&\vdots\\\
{{s_{n1}}}&\cdots&{{s_{nn}}}\end{pmatrix}\begin{pmatrix}{{A_{1}}}\\\ \vdots\\\
{{A_{n}}}\end{pmatrix}.$ (5)

If the transform matrix $S$ is invertible, the inverted transformation is

$\begin{pmatrix}{{A_{1}}}\\\ \vdots\\\
{{A_{n}}}\end{pmatrix}={S^{-1}}\begin{pmatrix}{{C_{1}}}\\\ \vdots\\\
{{C_{n}}}\end{pmatrix}.$ (6)

As the number $C_{i}$ is the scattering amplitude in channel $i$ and $A_{j}$
is the radiation amplitude from source $j$, all the scattering fields can be
controlled artificially according to Eq. 6.

There is also a remaining theoretical question to explain. What is the
reversible condition for the transform matrix? According to the linear algebra
theory, if $\\{\mathbf{h}_{1},\dots,\mathbf{h}_{n}\\}$ and
$\\{\mathbf{H}^{({\text{s}})}_{1},\dots,\mathbf{H}^{({\text{s}})}_{n}\\}$ are
respectively two set of linear independent functions, then $S$ is invertible.
First, the channel modes $\\{\mathbf{h}_{1},\dots,\mathbf{h}_{n}\\}$ are
solved from Helmholtz equations with different wavevectors, they of course are
linear independent each other. Second, we demonstrate with some examples that
$\\{\mathbf{H}^{({\text{s}})}_{1},\dots,\mathbf{H}^{({\text{s}})}_{n}\\}$ are
independent of each other when the sources are set correctly.

Example 1, two dipoles perpendicular to each other on metal surface as shown
is Fig. 1 (b) and (c), where $\mathbf{H}^{({\text{s}})}_{1}\neq
c\mathbf{H}^{({\text{s}})}_{2}$, $c$ is an arbitrary constant. In addition to
the above perpendicular electric dipoles, the mutually perpendicular magnetic
dipoles (or electric quadrupoles), etc., as the sources can excite independent
fields. In practice, although a scatterer illuminated by incident light
contains many sources, but it is difficult to independently control these
sources. Nevertheless, the same type of sources set at different positions can
make $\mathbf{H}^{({\text{s}})}_{1}\neq c\mathbf{H}^{({\text{s}})}_{2}$. The
following case is an example.

Example 2, two same electric dipoles $\mathbf{d}_{z1}$ and $\mathbf{d}_{z2}$
with different initial phase are set on the metal surface at different points,
just as sketched in Fig. 2. They also have the equation
$\mathbf{H}^{{\text{s}}}_{1}\neq c\mathbf{H}^{{\text{s}}}_{2}$.

Figure 2: The unidirectional scattering patterns excited by two vertical
dipoles in 2/2 structures. (a) for metal surface system. (b) for dielectric
slab system.

However, the key step to achieve the complete control in $n$-channel structure
is to set the $n$ sources correctly to excite $n$ independent fields. In this
case, the set of $n$ sources is referred to as a complete set of sources. In
other words, the condition of complete control is a complete set of sources.

In practice, the condition of $|S|\neq 0$ can guarantee the $n$ sources form a
complete set. If $|S|=0$, these sources should be rearranged. The mathematical
expression of the complete set of sources is

$\operatorname{rank}S=n.$ (7)

The above discussion is about the case of n sources set in n-channel
structure, which is referred to as $n/n$ structure. If $m$ sources are set in
a $n$-channel structure ($m/n$ structure), what should happen?

If $m>n$, Eq. 5 acts as $n$ equations, with $m$ elements $\\{A_{i}\\}$. $S$ is
a $n\times m$ matrix, if $\operatorname{rank}S=n$, the solutions $\\{A_{i}\\}$
always existed regardless of whatever values the $\\{C_{i}\\}$ take. The
complete control can be achieved. Here $\operatorname{rank}S=n$ means that,
there are $n$ sources out of m sources formed a complete set. So the complete
control in $m/n$ structure is also expressed with Eq. 7.

When $m<n$, We can select $m$ equations from Eq. 5, its coefficient matrix
$S_{m}$ is a $m\times m$ order submatrix. There are two cases are discussed
below.

If $|S_{m}|\neq 0$, $S_{m}$ is an invertible matrix, the corresponding $C_{i}$
(named selected $C_{i}$) contained in the selected equations can be assigned
with any values, and $A_{1},\dots,A_{m}$ are solved from these equations. The
other $C_{j}$ must be determined by the remaining equations of Eq. 5. It means
that, the waves in the $m$ channels can be controlled, this is the so-called
$m/n$ control, which needs $\operatorname{rank}S=m$.

If $|S_{m}|=0$, it can not guarantee the selected equations have solutions in
the case of that the selected $m$ $C_{i}$ take arbitrary values. So the $m/n$
control fails.

The both cases of $|S_{m}|=0$ and $|S_{m}|\neq 0$ may occur in a same $m/n$
structure even the rank $S$ equals to $m$, it depends on the selection of
equations. Therefore, in practice, one should first select the wanted channels
in which the waves are controllable, and then to set the $m$ sources correctly
to insure them independent. In this case the $m/n$ control can be achieved. So
in a $m/n$ structure, the most optimal control is $m/n$ control, If the source
setting is not appropriate, the $m/n$ control is not available.

Therefore, we can understand these works about four-channel structures[16, 17,
18, 19, 20]. As the circular polarized incident light contains only two
sources, only two channels can be selected as controllable channels. In other
words, the most optimal control should be 2/4 control, the 4/4 control is
unachievable.

Now we turn to show the applications of this theory in some structures. The
most typical PHE is sketched in Fig 1, but we would like to design another
novel and simple structure. Assume the same dipoles $\mathbf{d}_{z1}$ and
$\mathbf{d}_{z2}$ stand on metal surface at $x=-a$ and $x=a$, respectively
(2/2 structure). The channel 1($\mathbf{h}_{1}$) is in the region of $x<-a$
and channel 2 ($\mathbf{h}_{2}$) in $x>a$. If the dipole $\mathbf{d}_{z}$ is
set at $x=0$, $z=0$ the SPP wave can be written as
${\mathbf{H}^{(\text{s})}}=c\tfrac{x}{|x|}{d_{z}}{e^{i{k_{s}}\left|x\right|-\beta
z}}\mathbf{\hat{y}}$[17], here
$k_{s}=k\sqrt{\varepsilon_{\text{metal}}/(\varepsilon_{\text{metal}}+\varepsilon_{\text{air}})}$,
$k$ is the wavevector in vacuum (or air), $\beta=\sqrt{k_{s}^{2}-k^{2}}$, $c$
is a constant when the metal, the distance h of the dipole apart from the
metal surface, and incident wavelength $\lambda$ are determined. Therefore,
when taking $d_{z}=1$ and shifting it to $x=-a$, the $\mathbf{d}_{z1}$ should
excites the wave
${\mathbf{H}^{(\text{s})}_{1}}=c\tfrac{x+a}{|x+a|}{e^{i{k_{s}}\left|x+a\right|-\beta
z}}\mathbf{\hat{y}}$. In channel-1,
${\mathbf{H}^{(\text{s})}_{1}}=-c{e^{i{k_{s}}\left|x+a\right|-\beta
z}}\mathbf{\hat{y}}$, if we define $\mathbf{h}_{1}=-e^{-ik_{x}(x+a)-\beta
z}\mathbf{\hat{y}}$(amplitude is 1), the element
${s_{11}}=H_{1}^{(s)}/{h_{1}}=c$; in channel-2,
${\mathbf{H}^{(\text{s})}_{1}}=c{e^{i{k_{s}}(x+a)-\beta
z}}\mathbf{\hat{y}}=ce^{2ik_{s}a}{e^{i{k_{s}}(x-a)-\beta z}}\mathbf{\hat{y}}$,
$\mathbf{h}_{2}={e^{i{k_{s}}(x-a)-\beta z}}\mathbf{\hat{y}}$ so
${s_{21}}=c{e^{i2{k_{s}}a}}$. For the same reason, the dipole
$\mathbf{d}_{z2}$ contributes ${s_{12}}=c{e^{i2{k_{s}}a}}$, ${s_{22}}=c$. Then

${S^{-1}}=\frac{1}{{c(1-{e^{i2{k_{s}}a}})}}\begin{pmatrix}1&{-{e^{i{k_{s}}a}}}\\\
{-{e^{i{k_{s}}a}}}&1\end{pmatrix}.$ (8)

If we want the scattering wave to transmit to right, it means $C_{1}=0$;
$C_{2}=1$, and then Eqs. 6 and 8 gives

${A_{1}}=\frac{{-{e^{i{k_{s}}a}}}}{{c(1-{e^{i2{k_{s}}a}})}},{\text{
}}{A_{2}}=\frac{1}{{c(1-{e^{i2{k_{s}}a}})}}.$ (9)

As the above scattering wavefunctions $\mathbf{H}^{(\text{s})}_{1}$ and
$\mathbf{H}^{(\text{s})}_{2}$ are written in the case of $d_{z1}=d_{z2}=1$.
Equation 9 means $d_{z1}=A_{1}$ and $d_{z2}=A_{z2}$. In other words, the ratio
of ${d_{z2}}/{d_{z1}}={A_{2}}/{A_{1}}=-{e^{-i{k_{s}}a}}$ can cause the
scattering to right. For example, when assuming
${\varepsilon_{\text{metal}}}=-10+0.02i$, $h=0.2\lambda$ and $a=0.6\lambda$,
the SPP wavenumber is ${k_{s}}=(1.0541+0.0001171i)k$, then
${d_{2}}/{d_{1}}=0.09347-0.99474i$. The field distribution calculated under
this condition is displayed in Fig. 2(a), which is consistent with our
expectation.

By the way, we consider another 2/2 structure, the metal surface is replaced
by a dielectric slab with the thickness of $0.17\lambda$ and dielectric
constant of $9+0.02i$, $a=0.39\lambda$, $h=0.1\lambda$ and the waveguide
wavenumber is ${k_{s}}=(1.6108+0.003242i)k$. By using the same way, we
obtained ${d_{2}}/{d_{1}}=0.03954-0.98344i$. In this case, the unidirectional
scattering result is shown in Fig. 2(b).

A 4-channel structure can be formed by two crossed metal slits (Fig. 3(a)).
The four dipoles $\mathbf{d}_{1}$, $\mathbf{d}_{2}$, $\mathbf{d}_{3}$, and
$\mathbf{d}_{4}$ are set in different positions in different directions, as
shown in Fig. 3(e). This creates a 4/4 structure. It’s better to get S by
numerical simulation method. Since this structure has the $\pi/4$ rotational
symmetry, the simulation result from one dipole is sufficient. For instance,
excluding the dipoles $\mathbf{d}_{2}$, $\mathbf{d}_{3}$, and
$\mathbf{d}_{4}$, only $\mathbf{d}_{1}$ stands at the original position, the
excited SPP waves in channel-1 ($x>0$ direction ), channel-2 ($z>0$
direction), channel-3 ($x<0$ direction), and channel-4 ($z<0$ direction) are
denoted with $A$, $B$, $C$, $D$, respectively. Thus $s_{11}=A$, $s_{21}=B$,
$s_{31}=C$, and $s_{41}=D$. As the symmetry, the isolated $\mathbf{d}_{2}$
will contribute the result of $s_{12}=D$, $s_{22}=A$, $s_{32}=B$, and
$s_{42}=C$. And so on, the transform matrix is

$S={\begin{pmatrix}A&D&C&B\\\ B&A&D&C\\\ C&B&A&D\\\ D&C&B&A\end{pmatrix}},$
(10)

when assuming $\varepsilon_{\text{metal}}=-20+0.02i$ ,$a=0.1\lambda$ and the
slit width is $w=0.4\lambda$ . In these slits the SPP wavenumber is
${k_{s}}=(1.0965+0.0000539969i)k$. We can get $A=-0.58683-0.80182i$,
$B=D=0.02201+0.29451i$ and $C=-0.18300+0.83566i$ then $S$ and $S^{-1}$ are
obtained.

For this four-channel structure, there are four typical scattering
distribution patterns, they are shown in Fig. 3(a), (b), (c), and (d), and
named as a, b, c, and d patterns, respectively. For the a-pattern, $C_{1}=1$
and $C_{2}=C_{3}=C_{4}=0$, then $A_{1},\cdots,A_{4}$ are obtained from Eq. 6.
With the results of Eq. 6, the actual simulation result is shown in Fig. 3(a).
With the same way, the other results are displayed in Fig. 3(b), (c), and (d),
respectively.

Figure 3: The four typical scattering patterns in 4/4 structures. (a)-(d) for
the structure consisting of two crossed metal slits. The four dipoles are set
up as shown in (e). (f)-(i) for the structure composed of two crossed
dielectric slabs, the four dipoles are set up as shown in (j).

Another four-channel structure can be composed with two crossed dielectric
slabs, the four dipoles are set as shown in Fig. 3(j). Assuming the dielectric
constant is $9+0.02i$, $w=0.17\lambda$, and $a=0.1\lambda$. The a, b, c, d
patterns are obtained, and put in Fig. 3(f), (g),(h), and (i), respectively.

In summary, the scattering control theory in multi-channel structures are
proposed. The complete control requires a complete set of sources. The
application of this theory in a $n/n$ structure includes the following four
steps: first, determining the eigenmode in each channel; second, determining
the scattering waves excited by each source with unit or standard amplitude;
third, obtaining the transform matrix; and fourth, acquiring and setting the
amplitude of each source though the transform formula (Eq. 6) according to our
control intention. As the application demonstrations, the complete controls
for two-channel structures and four-channel structures are realized
theoretically. For a $m/n$ structure ($m<n$), one should first select the
wanted channels and set the sources correctly, then the $m/n$ control can be
realized. When $m>n$, the complete control needs at least $n$ sources to form
a complete set.

## References

  * Lee  _et al._ [2012] S.-Y. Lee, I.-M. Lee, J. Park, S. Oh, W. Lee, K.-Y. Kim, and B. Lee, Role of magnetic induction currents in nanoslit excitation of surface plasmon polaritons, Physical review letters 108, 213907 (2012).
  * Rodríguez-Fortuño  _et al._ [2013] F. J. Rodríguez-Fortuño, G. Marino, P. Ginzburg, D. O’Connor, A. Martínez, G. A. Wurtz, and A. V. Zayats, Near-field interference for the unidirectional excitation of electromagnetic guided modes, Science 340, 328 (2013).
  * Rodríguez-Fortuño  _et al._ [2014] F. J. Rodríguez-Fortuño, D. Puerto, A. Griol, L. Bellieres, J. Martí, and A. Martínez, Sorting linearly polarized photons with a single scatterer, Optics letters 39, 1394 (2014).
  * Rodriguez-Fortuno  _et al._ [2014] F. J. Rodriguez-Fortuno, I. Barber-Sanz, D. Puerto, A. Griol, and A. Martínez, Resolving light handedness with an on-chip silicon microdisk, ACS Photonics 1, 762 (2014).
  * Mitsch  _et al._ [2014] R. Mitsch, C. Sayrin, B. Albrecht, P. Schneeweiss, and A. Rauschenbeutel, Quantum state-controlled directional spontaneous emission of photons into a nanophotonic waveguide, Nature communications 5, 1 (2014).
  * Petersen  _et al._ [2014] J. Petersen, J. Volz, and A. Rauschenbeutel, Chiral nanophotonic waveguide interface based on spin-orbit interaction of light, Science 346, 67 (2014).
  * Kapitanova  _et al._ [2014] P. V. Kapitanova, P. Ginzburg, F. J. Rodríguez-Fortuño, D. S. Filonov, P. M. Voroshilov, P. A. Belov, A. N. Poddubny, Y. S. Kivshar, G. A. Wurtz, and A. V. Zayats, Photonic spin hall effect in hyperbolic metamaterials for polarization-controlled routing of subwavelength modes, Nature communications 5, 1 (2014).
  * Espinosa-Soria and Martinez [2016] A. Espinosa-Soria and A. Martinez, Transverse spin and spin-orbit coupling in silicon waveguides, IEEE Photonics Technology Letters 28, 1561 (2016).
  * Coles  _et al._ [2016] R. Coles, D. Price, J. Dixon, B. Royall, E. Clarke, P. Kok, M. Skolnick, A. Fox, and M. Makhonin, Chirality of nanophotonic waveguide with embedded quantum emitter for unidirectional spin transfer, Nature communications 7, 1 (2016).
  * Picardi  _et al._ [2017] M. F. Picardi, A. Manjavacas, A. V. Zayats, and F. J. Rodríguez-Fortuño, Unidirectional evanescent-wave coupling from circularly polarized electric and magnetic dipoles: An angular spectrum approach, Phys. Rev. B 95, 245416 (2017).
  * Bliokh  _et al._ [2014] K. Y. Bliokh, A. Y. Bekshaev, and F. Nori, Extraordinary momentum and spin in evanescent waves, Nature communications 5, 1 (2014).
  * Bliokh and Nori [2015] K. Y. Bliokh and F. Nori, Transverse and longitudinal angular momenta of light, Physics Reports 592, 1 (2015).
  * Bliokh  _et al._ [2015a] K. Y. Bliokh, F. J. Rodríguez-Fortuño, F. Nori, and A. V. Zayats, Spin–orbit interactions of light, Nature Photonics 9, 796 (2015a).
  * Bliokh  _et al._ [2015b] K. Y. Bliokh, D. Smirnova, and F. Nori, Quantum spin hall effect of light, Science 348, 1448 (2015b).
  * Van Mechelen and Jacob [2016] T. Van Mechelen and Z. Jacob, Universal spin-momentum locking of evanescent waves, Optica 3, 118 (2016).
  * Picardi  _et al._ [2018] M. F. Picardi, A. V. Zayats, and F. J. Rodríguez-Fortuño, Janus and huygens dipoles: Near-field directionality beyond spin-momentum locking, Phys. Rev. Lett. 120, 117402 (2018).
  * Gao  _et al._ [2021] Q. Gao, Y.-S. Zhou, and L.-M. Zhao, Continuous photon hall effect on metal surface, EPL (Europhysics Letters) 133, 27001 (2021).
  * Luxmoore  _et al._ [2013] I. J. Luxmoore, N. A. Wasley, A. J. Ramsay, A. C. T. Thijssen, R. Oulton, M. Hugues, S. Kasture, V. G. Achanta, A. M. Fox, and M. S. Skolnick, Interfacing spins in an ingaas quantum dot to a semiconductor waveguide circuit using emitted photons, Phys. Rev. Lett. 110, 037402 (2013).
  * Wang  _et al._ [2017] Y.-H. Wang, R.-C. Jin, J.-Q. Li, F. Zhong, H. Liu, I. Kim, Y. Jo, J. Rho, and Z.-G. Dong, Photonic spin hall effect by the spin-orbit interaction in a metasurface with elliptical nano-structures, Applied Physics Letters 110, 101908 (2017).
  * Zhao and Zhou [2021] L.-M. Zhao and Y.-S. Zhou, Tunable multichannel photonic spin hall effect in metal-dielectric-metal waveguide, Scientific Reports 11, 1 (2021).


#  A Survey Selection Correction using Nonrandom Followup with an Application
to the Gender Entrepreneurship Gap

Clint Harris111 University of Wisconsin-Madison<EMAIL_ADDRESS>Jon
Eckhardt222 University of Wisconsin-Madison<EMAIL_ADDRESS>Brent
Goldfarb333 University of Maryland<EMAIL_ADDRESS>
###### Abstract

Selection into samples undermines efforts to describe populations and to
estimate relationships between variables. We develop a simple method for
correcting for sample selection that explains differences in survey responses
between early and late respondents with correlation between potential
responses and preference for survey response. Our method relies on researchers
observing the number of data collection attempts prior to each individual’s
survey response rather than covariates that affect response rates without
affecting potential responses. Applying our method to a survey of
entrepreneurial aspirations among undergraduates at University of Wisconsin-
Madison, we find suggestive evidence that the entrepreneurial aspiration rate
is larger among survey respondents than the population, as well as the male-
female gender gap in the entrepreneurial aspiration rate, which we estimate as
21 percentage points in the sample and 19 percentage points in the population.
Our results suggest that the male-female gap in entrepreneurial aspirations
arises prior to direct exposure to the labor market.

JEL Codes: C83, L26, J16
Keywords: Survey Nonresponse, Attrition, Sample Selection, Entrepreneurship,
Gender

##  1 Introduction

Estimates of relationships between variables in a sample of a population are
informative about the relationships between those variables in the population
of interest when the sample is a random draw from the population (Fisher,
1922). To take advantage of this remarkable property, researchers commonly
choose whether an individual in a population will be surveyed at random.
Unfortunately, surveyed individuals choose whether they will respond to
surveys nonrandomly. This nonrandom choice breaks the connection between
between in-sample estimates and population estimands by introducing selection
bias, which undermines inference about related research questions (Heckman,
1979).

This paper introduces a method for correcting for survey nonresponse by
extrapolating from differences between early and late survey respondents in
measured variables to differences between late and never responders in those
variables. It has modest data demands, not requiring observation of baseline
characteristics that predict survey nonresponse, unlike covariate adjustment
based methods such as those described by Molina Millán and Macours
(2017).444The method of this paper easily incorporates covariate adjustment if
variables correlated with survey response are observed. It also places minimal
demands on the survey administrator, not requiring randomization of incentives
or followup intensity, as in methods described by Dutz, Huitfeldt, Lacouture,
Mogstad, Torgovitsky, and Van Dijk (2021) and DiNardo, Matsudaira, McCrary,
and Sanbonmatsu (2021), respectively.555Ex post randomization of which late
responders’ surveys to discard in a setting with multiple reminders for
everyone would present a scenario in which the method of DiNardo, Matsudaira,
McCrary, and Sanbonmatsu (2021) is directly applicable, with a slight tweak in
the interpretation of who intensive followup compliers are. Our method relies
on the same intuition without reductions in sample size from randomly
discarding data (or never receiving it from followup compliers who are
randomized into receiving low followup intensity).  The present method may
therefore be particularly attractive in applications using secondary data,
where a researcher’s estimation goals did not inform the data collection.
Furthermore, it does not require financial or similar incentives, reducing
costs (especially in large surveys) and avoiding the possibility of incentives
crowding out intrinsic motivation for some individuals, causing them to reduce
their response effort rather than increase it (Gneezy and Rustichini, 2000).

Our method shares all of the above advantages with that of Behaghel, Crépon,
Gurgand, and Le Barbanchon (2015) for applications focused on treatment effect
estimation. The modeling distinction between these methods is that ours
requires assumptions (e.g. normality) on the joint distributions of unobserved
determinants of the outcome of interest and survey response aversion that are
allowed to vary by treatment status, while theirs instead assumes that the
distribution of survey response aversion for treated individuals is a
monotonic transformation of the distribution for untreated individuals. Our
alternative distributional assumptions lead to three key distinctions between
results produced by our method and those produced by the method of Behagel et
al. 2015. First, our method provides population level conditional expectation
estimates, so it estimates average treatment effects for the population in
treatment effect applications, rather than local average treatment effects
specific to individuals with low survey response aversion. This distinction
also extends our method’s applicability to estimation of descriptive
quantities such as population averages and comparisons between nonrandomly
assigned groups (e.g. gender gaps), which may not have policy-relevant local
averages. Second, our method produces point estimates of quantities of
interest rather than bounds. Third, and finally, our method allows for
heterogeneity between groups (such as by treatment status) in the correlation
between outcomes and survey response aversion, which can drive differential
selection bias between groups even in settings where they have similar survey
response rates.666Our method is robust to settings in which a treatment with
no effect on outcomes increases (decreases) survey responses for high outcome
(low outcome) individuals, which would likely yield positive effect estimates
without correction or with a correction that relied on monotonic relationships
between survey responses and outcomes between groups.

We use our method to estimate the gender-gap in the prevalence of
entrepreneurial aspirations among college students, using a survey of the
undergraduate population of University of Wisconsin-Madison. A large body of
work documents substantial differences in rates of entrepreneurship between
men and women (Guzman and Kacperczyk, 2019), and suggests that women’s
underrepresentation in entrepreneurship may harm consumers if valuable
products and services are overlooked by male entrepreneurs. Explanations for
women’s underrepresentation in entrepreneurship can be broadly categorized
into nature vs. nurture, with nature-based explanations focusing on gender
differences in relevant intrinsic preferences, such as for riskand
competition, and with nurture-based explanations focusing on gender
differences in external factors, such as labor market discrimination,
differences in social norms surrounding childcare, and premarket disparities
in education or encouragement relevant to entrepreneurship. Our nonresponse-
corrected estimate of the entrepreneurial aspirations gender gap is large (19
percentage points) and of similar magnitude to the uncorrected estimate (21
percentage points), suggesting that premarket factors are major determinants
of the entrepreneurship gender gap.

##  2 Survey Nonresponse Correction with Multiple Reminders

We consider a randomly surveyed sample of $N$ individuals indexed by $i$ who
receive $T$ requests to fill out a survey, with survey requests indexed by
$t$. We denote individual $i$’s potentially unobserved outcome of interest as
$Y_{i}$ and the binary indicator for them responding to request $t$ or a prior
request as $S_{it}$, such that $S_{iT}=1$ if they ever respond and $S_{iT}=0$
if they never respond. We are interested in estimation of quantities such as
the unconditional expectation of the outcome in the population,
$\mathbb{E}[Y_{i}]$, conditional expectations of the outcome for various
groups defined by values of observed covariates such as
$\mathbb{E}[Y_{i}|W_{i}=w]$, and differences between groups in their expected
outcomes such as
$\mathbb{E}[Y_{i}|W_{i}=w]-\mathbb{E}[Y_{i}|W_{i}=w^{\prime}]$, where
$W_{i}\subseteq X_{i}$ and $X_{i}$ is the set of covariates observed for
surveyed individuals.777The comparison
$\mathbb{E}[Y_{i}|W_{i}=w]-\mathbb{E}[Y_{i}|W_{i}=w^{\prime}]$ gives the
average treatment effect of $W_{i}=w$ versus $W_{i}=w^{\prime}$ in cases in
which $W_{i}$ is randomized, with similar quantities being relevant to other
research designs focused on measuring treatment effects. The application in
this paper investigating the gender gap in entrepreneurial aspirations
estimates such a quantity, following procedures that would consistently
estimate the treatment effect of gender if it were randomly assigned.

###  2.1 Problem

In settings with survey nonresponse, the conditional sample mean among
individuals for whom $W_{i}=w$ is a consistent estimate of
$\mathbb{E}[Y_{i}|W_{i}=w,S_{iT}=1]$, which may differ from
$\mathbb{E}[Y_{i}|W_{i}=w]$ as described by Heckman (1979). The expected
nonresponse bias for group $W_{i}=w$ is given by

$B_{w}=\mathbb{E}[Y_{i}|W_{i}=w,S_{iT}=1]-\mathbb{E}[Y_{i}|W_{i}=w].$ (1)

The difference in sample means between groups defined by $W_{i}=w$ and
$W_{i}=w^{\prime}$ provides an estimate of

$\begin{split}&\mathbb{E}[Y_{i}|W_{i}=w,S_{iT}=1]-\mathbb{E}[Y_{i}|W_{i}=w^{\prime},S_{iT}=1]\\\
=&\mathbb{E}[Y_{i}|W_{i}=w]-\mathbb{E}[Y_{i}|W_{i}=w^{\prime}]+B_{w}-B_{w^{\prime}},\end{split}$
(2)

which is only equal to the population difference between groups if their
survey nonresponse bias terms are equal. In settings with 100% survey response
rates (or random nonresponse), each of these bias terms is equal to zero. It
follows that in such cases we can consistently estimate
$\mathbb{E}[Y_{i}|W_{i}=w]$ for any $W_{i}\subseteq X_{i}$ and $w$ in the
support of $W_{i}$ by calculating the relevant conditional sample mean, and we
can compare or average these quantities among groups as desired.

###  2.2 Solution

It is hard to justify a priori assumptions about nonresponse bias magnitudes
or random nonresponse, so we propose a method for extrapolating from observed
values of $Y_{i}$ for individuals who respond to the survey to the rest of the
population. In our application, outcomes of interest are binary choices made
by agents, so we assume that

$Y_{i}=\begin{cases}&1\mbox{ if }X_{i}\beta+\epsilon_{i}\geq 0,\\\ &0\mbox{
otherwise,}\end{cases}$ (3)

wherein $X_{i}$ is a $1\times K$ vector containing a constant and observed
variables that predict outcomes, $\beta$ is an $K\times 1$ vector of marginal
effects of variables in $X_{i}$, and $\epsilon_{i}$ is a mean zero
residual.888The case in which $Y_{i}$ is continuous will be described in a
future Appendix. Importantly, the outcome process in (3) contains no $t$
subscripts, so effects of notifications (or time between them) on the outcome
or its measurement at various points in time are assumed to be zero. We
similarly assume that individuals make a binary choice regarding responding to
survey request $t$ or a prior request according to the rule

$S_{it}=\begin{cases}&1\mbox{ if }Z_{i}\alpha_{t}+u_{i}\geq 0,\\\ &0\mbox{
otherwise,}\end{cases}$ (4)

wherein $Z_{i}$ is a $1\times L$ vector containing a constant and observed
variables that predict outcomes, $\alpha_{t}$ is an $L\times 1$ vector of
marginal effects of variables in $Z_{i}$, and $u_{i}$ is a mean zero
residual.999A similar model would treat the response decision as ordered and
would assume individual $i$ responds after notification $t$ and before
notification $t+1$ if $\mu^{t}\leq Z_{i}\alpha+u_{i}<\mu^{t+1}$ with
normalizations $\mu^{1}=0$ and $\mu^{T+1}=\inf$. The model in (4) is
attractive because it allows for individuals with different values of $Z_{i}$
to differ in their intensities of response to reminders, essentially using
$\mu^{t}(Z_{i})=Z_{i}\mu^{t}$ in the ordered model and defining
$\alpha_{t}=\alpha-\mu^{t}$.

We follow Heckman (1979) and Van de Ven and Van Praag (1981), who extended
Heckman’s selection correction method to the case with a binary outcome, in
assuming that unobserved residuals of binary outcome equations are normally
distributed, in our case as

$\begin{bmatrix}\epsilon_{i}\\\
u_{i}\end{bmatrix}\sim\mathcal{N}\left(\begin{bmatrix}0\\\
0\end{bmatrix},\begin{bmatrix}1\qquad\rho\\\ \rho\qquad
1\end{bmatrix}\right).$ (5)

We discuss the implications of these assumptions relative to some alternatives
in Section 2.3, and we describe tests of these assumption that can be
performed in applications with more than one reminder in Section 2.4. Under
this assumption, the log-likelihood for the sample is

$\begin{gathered}\mathcal{L}\left(\beta,\alpha_{1},...,\alpha_{T},\rho|X_{i},Z_{i}\right)=\\\
\sum_{t=1}^{T}\sum_{i=1}^{N}S_{it}(Y_{i}\ln(\Phi_{2}(X_{i}\beta,Z_{i}\alpha_{t},\rho))+(1-Y_{i})\ln(\Phi_{2}(-X_{i}\beta,Z_{i}\alpha_{t},-\rho)))\\\
+(1-S_{it})\ln(1-\Phi(Z_{i}\alpha_{t})),\end{gathered}$ (6)

wherein $\Phi(\cdot)$ is the standard normal CDF and $\Phi_{2}(\cdot)$ is the
standard bivariate normal CDF. The values
$(\hat{\beta},\hat{\alpha}_{1},...,\hat{\alpha}_{T},\hat{\rho})$ that maximize
this likelihood are estimates of their corresponding population parameters,
and the estimated conditional probability of $Y_{i}$ given $X_{i}$ is
$\Phi(X_{i}\hat{\beta})$.

It is helpful to think through the estimator’s performance in a few extreme
cases to build intuition. First, in the case where $t=1$ and $X_{i}=Z_{i}$,
maximizing the likelihood subject to the constraint $\rho=0$ provides the same
estimates of $\beta$ as a probit regression of $Y_{i}$ on $X_{i}$, which is
valid if survey nonresponse is random conditional on $X_{i}$. Without an
unconstrained $\rho$ parameter, the estimator reduces to the probit version of
Heckman’s sample selection correction method, which is reliable for $T=1$ in
cases where there is at least one variable in $Z_{i}$ that is not also in
$X_{i}$. This motivates the use of randomized survey response incentives or
followup intensity as advocated by Dutz, Huitfeldt, Lacouture, Mogstad,
Torgovitsky, and Van Dijk (2021) and DiNardo, Matsudaira, McCrary, and
Sanbonmatsu (2021). In cases where $t>1$ and $X_{i}=Z_{i}$, the model is
identified by differences in average values of $Y_{i}$ for individuals who
respond to the survey at different times. This is equivalent to Heckman’s
correction with $X_{i}\subset Z_{i}$ under an equivalent model that defines
$Z_{i}=[X_{i},X_{i}S_{i2},...X_{i}S_{iT}]$, with an $1\times KT$ coefficient
vector $\alpha$ present in survey response equations for all survey
notifications.

###  2.3 Assumption Violation Implications

It is useful to think through implications of violations of model assumptions.
First, if the outcome (or nonclassical error in its measurement) changes in
response to the passage of time or in response to survey reminders, at odds
with the time invariance assumption implied by (3), then the estimated model’s
attribution of differences between early and late responders to survey
response preferences will be inaccurate. Concerns relating to the passage of
time can be eliminated if survey responses are restricted to time invariant
quantities, which survey administrators can encourage by careful wording of
questions.101010To use a salient example, we would expect our method to
perform poorly if used during a pandemic on a survey asking individuals if
they currently have the illness because average values of this quantity would
change over time. We would expect it to perform better if used on a similar
survey asking individuals if they had the illness as of a specified date (such
as the date the survey was sent out). Direct effects of reminders on time-
invariant outcomes seem unlikely, but effects on their measurement could occur
if reminders affect survey effort. For instance, sending a large number of
reminders in quick succession could annoy respondents and cause them to
provide responses that are uncorrelated or even negatively correlated with
true values.

Second, if the distributional assumption on $\epsilon_{i}$ and $u_{i}$ in (5)
is inaccurate, the out of sample extrapolation will suffer. In cases with a
single reminder, estimates of conditional in-sample means given by
$\Phi(X_{i}\hat{\beta}+\hat{\rho}\lambda(Z_{i}\hat{\alpha}))$ (where
$\lambda(Z_{i}\alpha)$ is the inverse-Mills ratio evaluated at $Z_{i}\alpha$)
for early and late responders will be equal to sample means for both subgroups
(the constant and the parameter $\rho$ are chosen to ensure this), but out-of-
sample fit may be poor. In cases with more than one reminder, the in-sample
fit could suffer as well, which motivates the test of this assumption
described in section 2.4.

###  2.4 Testing

In settings with more than one reminder, the model is overidentified. Assuming
that the passage of time and the receipt of notifications are excludeable from
the outcome equation, the bivariate normality assumption on model residuals
given in (5) can be assessed by estimating the log-Likelihood

$\begin{gathered}\mathcal{L}\left(\beta,\alpha_{1},...,\alpha_{T},\rho_{2},...,\rho_{T}|X_{i},Z_{i}\right)=\sum_{t=1}^{T}\sum_{i=1}^{N}(1-S_{it})\ln(1-\Phi(Z_{i}\alpha_{t}))\\\
+\sum_{t=1}^{2}\sum_{i=1}^{N}S_{it}\left(Y_{i}\ln(\Phi_{2}(X_{i}\beta,Z_{i}\alpha_{t},\rho_{2}))+(1-Y_{i})\ln(\Phi_{2}(-X_{i}\beta,Z_{i}\alpha_{t},-\rho_{2}))\right)\\\
+\sum_{t=3}^{T}\sum_{i=1}^{N}S_{it}\left(Y_{i}\ln(\Phi_{2}(X_{i}\beta,Z_{i}\alpha_{t},\rho_{t}))+(1-Y_{i})\ln(\Phi_{2}(-X_{i}\beta,Z_{i}\alpha_{t},-\rho_{t}))\right),\\\
\end{gathered}$ (7)

and testing the hypothesis $\hat{\rho}_{2}=\hat{\rho}_{t}$ for $t=3,...,T$.
This test essentially compares the correlation parameter, $\hat{\rho}_{2}$,
that rationalizes the differences in outcomes between individuals who respond
to the first or second notification to the correlation parameters that
rationalize outcomes for those who respond to notification $t$ or before, for
all $t>2$, given the estimate of $\beta$ from the first two
notifications.111111For instance, when $X_{i}=Z_{i}=1$, $\alpha_{t}$ is
exactly identified by the survey response rate up to reminder $t$, $\beta$ is
identified as the predicted population expectation of $Y_{i}$ using only the
first two notifications, and $\rho_{t}$ for $t>2$ rationalizes the within-
sample average of $Y_{i}$ up to notification $t$ taking $\hat{\beta}$ as given
from the first two notifications. The equality of these correlation parameters
is necessary, but not sufficient, to determine that the true distribution of
errors is bivariate normal.

This test is unduly strict for some applications. It will reject the null of
bivariate normality in cases in which the correlation parameter that
rationalizes outcomes at $t=1$ and $t=2$ fails to rationalize outcomes at some
larger values of $t$, even if it does accurately rationalize later outcomes on
average. If the model’s prediction errors for outcomes from estimates only
using the first two notifications are independent of $t$ in and out-of-sample,
model estimates of population expectations will be consistent, with
reliability increasing in the number of reminders (such that the $\hat{\rho}$
of best fit averages across the prediction errors rather than fitting
them).121212This strictness is similar to that of pre-trend tests used in
difference-in-differences applications that check whether all pre-treatment
deviations from parallel trends are equal to zero. In such cases estimates of
policy effects that average over a large number of post-treatment time periods
are valid if parallel trends violations are independent of time.  A more
lenient test can be implemented by maximizing

$\begin{gathered}\mathcal{L}\left(\beta,\alpha_{1},...,\alpha_{T},\rho_{2},...,\rho_{T}|X_{i},Z_{i}\right)=\sum_{t=1}^{T}\sum_{i=1}^{N}(1-S_{it})\ln(1-\Phi(Z_{i}\alpha_{t}))\\\
+\sum_{t=1}^{2}\sum_{i=1}^{N}S_{it}\left(Y_{i}\ln(\Phi_{2}(X_{i}\beta,Z_{i}\alpha_{t},\rho_{2}))+(1-Y_{i})\ln(\Phi_{2}(-X_{i}\beta,Z_{i}\alpha_{t},-\rho_{2}))\right)\\\
+\sum_{i=1}^{N}S_{iT}\left(Y_{i}\ln(\Phi_{2}(X_{i}\beta,Z_{i}\alpha_{T},\rho_{T}))+(1-Y_{i})\ln(\Phi_{2}(-X_{i}\beta,Z_{i}\alpha_{T},-\rho_{T}))\right),\\\
\end{gathered}$ (8)

and testing the hypothesis $\hat{\rho}_{2}=\hat{\rho}_{T}$. This test checks
whether the estimated conditional expectations of $Y_{i}$ and correlation
parameter obtained using only the first two notifications also rationalizes
average observed outcomes among all individuals who respond to the survey,
without penalizing prediction error for distinct response groups.

##  3 Application: Gender Entrepreneurship Gap

We use our method to estimate the gender gap in entrepreneurial aspirations
among college students. There is a large gender gap in entrepreneurship
between working age adults (Aldrich, 2005), driven at least in part by a
gender gap in early stage funding for new ventures (Canning, Haque, and Wang,
2012; Greene, Hart, Gatewood, Brush, and Carter, 2003). If there is no such
gap in intentions among individuals prior to entering the labor market, we can
conclude that the gap is due to market forces and not gender differences in
other factors such as personal preferences relevant to entrepreneurship or
premarket discrimination. We contribute to this literature by estimating the
entrepreneurial aspiration gender gap using the method described in the
previous section to address selection bias.

###  3.1 Data

We use a survey of the undergraduate population of University of Wisconsin -
Madison that was implemented every Fall from 2015 to 2022, as well as in
Spring 2020 and 2021. We match the survey data to administrative data on
students, which contains information on individuals regardless of whether they
responded to the survey. This survey asked students whether they intended to
pursue a career in entrepreneurship. All responses were given as “Yes”, “No”,
or “I Don’t Know”, where we code responses as binary based on whether the
answer given was “Yes”. Summary variables for those surveyed broken down by if
and when they responded are provided in Table 1.

Table 1: Summary Statistics

| Early Respondents | Late Respondents | Nonrespondents
---|---|---|---
| Men | Women | Men | Women | Men | Women
| (1) | (2) | (3) | (4) | (5) | (6)
Intention | 0.374 | 0.167 | 0.363 | 0.179 | . | .
Business Major | 0.128 | 0.093 | 0.130 | 0.098 | 0.125 | 0.086
STEM Major | 0.455 | 0.276 | 0.455 | 0.269 | 0.432 | 0.261
ACT Math | 30.220 | 28.025 | 30.252 | 27.912 | 29.793 | 27.520
ACT Verbal | 55.718 | 56.478 | 54.747 | 55.667 | 54.239 | 54.798
Racial Minority | 0.242 | 0.227 | 0.287 | 0.264 | 0.302 | 0.281
International | 0.098 | 0.063 | 0.108 | 0.070 | 0.103 | 0.071
First-Gen College | 0.142 | 0.177 | 0.144 | 0.180 | 0.154 | 0.193
In-State | 0.617 | 0.649 | 0.576 | 0.609 | 0.546 | 0.562
Year 2 Student | 0.220 | 0.219 | 0.229 | 0.230 | 0.248 | 0.258
Year 3 Student | 0.176 | 0.183 | 0.184 | 0.191 | 0.236 | 0.245
Year 4 Student | 0.168 | 0.181 | 0.154 | 0.164 | 0.204 | 0.211
Year 5+ Student | 0.059 | 0.037 | 0.055 | 0.037 | 0.069 | 0.042
2016, Fall | 0.128 | 0.135 | 0.092 | 0.101 | 0.103 | 0.096
2017, Fall | 0.152 | 0.167 | 0.074 | 0.083 | 0.102 | 0.096
2018, Fall | 0.076 | 0.081 | 0.112 | 0.130 | 0.100 | 0.100
2019, Fall | 0.082 | 0.088 | 0.079 | 0.085 | 0.105 | 0.107
2020, Spring | 0.065 | 0.073 | 0.112 | 0.114 | 0.095 | 0.098
2020, Fall | 0.104 | 0.091 | 0.161 | 0.161 | 0.099 | 0.101
2021, Spring | 0.068 | 0.065 | 0.083 | 0.079 | 0.100 | 0.103
2021, Fall | 0.085 | 0.072 | 0.143 | 0.119 | 0.099 | 0.103
2022, Fall | 0.044 | 0.037 | 0.100 | 0.079 | 0.100 | 0.102
Observations | 9563 | 12407 | 8703 | 10152 | 106582 | 114979

_Notes:_ Variable means for men and women by responder group. Early
Respondents are individuals who respond prior to any reminders, while Late
Respondents are individuals who respond any time after the first reminder.

###  3.2 Estimated Quantities

We begin by estimating the rate of entrepreneurial intention among men and
women separately in the highest response rate term (Fall, 2020), without using
the predictor variables listed in Table 1. This parsimonious specification is
pedagogically attractive because it restricts attention to our methodological
contribution, while also clarifying the close connection of our method to a
Heckman correction, especially in the specific case of randomized followup
intensity as proposed by DiNardo, Matsudaira, McCrary, and Sanbonmatsu (2021).
Specifically, we maximize the log-Likelihood in (6) and calculate
$\Phi(\hat{\beta}_{0})$ separately for men and women, where $\hat{\beta}_{0}$
is the coefficient on the constant. We also calculate
$\Phi(\hat{\beta}_{0}+\hat{\rho}\Phi^{-1}(1-e_{i}))$ for the range of values
of $e_{i}\in[0,1]$ separately for men and women, in order to trace out
predicted entrepreneurial aspirations for individuals across percentiles of
survey response aversion. Importantly, we always estimate gender-specific
residual correlation ($\rho$) parameters.131313Allowing group specific
residual correlations is potentially important in applications involving
treatment effect estimation. For instance, in a randomized controlled trial
with post-experiment followup attrition, allowing $\rho$ to be different for
treated and untreated individuals can correct for the possibility of
differential selection bias, which is arguably the main threat to
identification in such settings.

We also estimate the average gender gap, and we perform an Oaxaca-Blinder
decomposition to identify variables that are relevant to the gender gap
(Oaxaca, 1973; Blinder, 1973). This method considers, for each covariate in
$X_{i}$ that determines outcomes, the extent to which the gender gap would be
reduced if women were like men with respect to that single variable, and how
much it would be reduced if women were like men with respect to how much that
variable determines outcomes (given by the element of $\beta$ corresponding to
the variable in a model where $X_{i}$).

The counterfactual in which women are “like men” with respect to a given
variable or its coefficient can be interpreted multiple ways. Specifically, we
define $F_{i}$ as a binary indicator for female gender, and we estimate
versions of outcome process (3) of the form

$Y_{i}=X_{i}\beta+(F_{i}\times X_{i})\delta+\epsilon_{i}$ (9)

where $\times$ defines the dot product such that $F_{i}\times X_{i}$ gives the
interaction of the female indicator and other variables, with $X_{i}$ here
defined differently from above as the variables other than gender that are
assumed to determine the outcome. For a single variable $W_{i}\subseteq
X_{i}$, we estimate the contribution to the average gender gap of that
variable as

$\hat{\Delta}_{W}=\sum_{i=1}^{N}\frac{\Phi(X_{i}\hat{\beta}+(F_{i}\times
X_{i})\hat{\delta}+F_{i}\times(\bar{W}^{M}-\bar{W}^{F})(\hat{\beta}_{W}+\hat{\delta}_{W}))-\Phi(X_{i}\hat{\beta}+(F_{i}\times
X_{i})\hat{\delta})}{N},$ (10)

wherein $\bar{W}^{M}$ is the sample mean of $W_{i}$ for men, $\bar{W}^{F}$ is
the sample mean of $W_{i}$ for women, $\hat{\beta}_{W}$ is the coefficient on
$W_{i}$ and $\hat{\delta}_{W}$ is the coefficient on the interaction between
$W_{i}$ and $F_{i}$.141414This is simpler than the commonly used extension of
the Oaxaca-Blinder Decomposition developed by Fairlie (2005) for nonlinear
models. Our method provides predictions (under conditional independence and
homogeneous effects assumptions on $W_{i}$) for the gap in $Y_{i}$ under
counterfactual policies that change $W_{i}$ the same amount for all women.
Fairlie’s method provides predictions (under the same assumptions) for the gap
in $Y_{i}$ under more nuanced counterfactual policies that shift $W_{i}$
different amounts for different women to match men’s marginal distribution of
$W_{i}$.  We similarly estimate the contribution to the average gender gap of
that variable’s differential effect on $Y_{i}$ between men and women as

$\hat{\Delta}_{\beta_{W}}=\sum_{i=1}^{N}\frac{\Phi(X_{i}\hat{\beta}+(F_{i}\times
X_{i})\hat{\delta}-(F_{i}\times
W_{i}^{F})\hat{\delta}_{W})-\Phi(X_{i}\hat{\beta}+(F_{i}\times
X_{i})\hat{\delta})}{N}.$ (11)

###  3.3 Results

We begin by presenting estimates of average rates of entrepreneurial intention
among men and women for Fall, 2020 without using the predictor variables
listed in Table 1. We present these results graphically in Figures 1 for men
and Figure 2 for women. The top panel of both figures shows means within each
response group, with individuals responding to the third or later reminder
pooled together for the graph, but not for the estimation. The top panel also
includes uncorrected and corrected extrapolations using within-sample means
(equivalent to results from an in-sample probit with only a constant) and our
method, respectively. The bottom panel of each figure adds extrapolated
estimates from the overidentification test from maximizing the likelihood in
(7), with overlapping confidence intervals for all extrapolations showing
that, within this one year, we are unable to reject the uncorrected method
null of no selection bias or the corrected method null of bivariate normality
on residuals for both men and women.

Figure 1: Entrepreneurial Intentions for Men

_Notes:_ Estimated average entrepreneurial aspiration rates. Uncorrected
estimates are calculated as the sample average among all survey respondents.
Corrected estimates are calculated as
$\Phi(\hat{\beta}_{0}+\hat{\rho}\Phi^{-1}(1-e_{i}))$ for the range of values
of survey aversion percentiles, $e_{i}\in[0,1]$.. Corrected extrapolations 1-2
and 1-3 provide estimates of
$\Phi(\hat{\beta}_{0}+\hat{\rho}_{2}\Phi^{-1}(1-e_{i}))$ and
$\Phi(\hat{\beta}_{0}+\hat{\rho}_{3}\Phi^{-1}(1-e_{i}))$, respectively,
estimated by maximizing the likelihood in (7).

Figure 2: Entrepreneurial Intentions for Women

_Notes:_ Estimated average entrepreneurial aspiration rates. Uncorrected
estimates are calculated as the sample average among all survey respondents.
Corrected estimates are calculated as
$\Phi(\hat{\beta}_{0}+\hat{\rho}\Phi^{-1}(1-e_{i}))$ for the range of values
of survey aversion percentiles, $e_{i}\in[0,1]$. Corrected extrapolations 1-2
and 1-3 provide estimates of
$\Phi(\hat{\beta}_{0}+\hat{\rho}_{2}\Phi^{-1}(1-e_{i}))$ and
$\Phi(\hat{\beta}_{0}+\hat{\rho}_{3}\Phi^{-1}(1-e_{i}))$, respectively,
estimated by maximizing the likelihood in (7).

We present estimates of the average entrepreneurial aspiration rate for men
and women for the entire population, using all years of our survey, in Table
2. This table also includes the estimated gender gaps in aspirations from the
uncorrected and corrected method, along with a decomposition of the
contributions of variables and their coefficients to the gender gap,
calculated via equations (10) and (11), respectively.

The panel denoted $X$ shows mean values for each variable in $X_{i}$ for both
men and women who responded to the survey in columns (1) and (2). Columns (4)
and (5) provide the sample means of each variable for men and women for all
individuals in our data regardless of whether they responded to the survey.
Columns (3) and (6) report estimated contributions of each variable in $X_{i}$
to the gender gap calculated as described in equation (10) using values of
$X_{i}$ for all individuals in our data, with the uncorrected estimates using
$\hat{\beta}$ from a within-survey probit and the corrected estimates coming
from maximizing the log-Likelihood in (6).

The panel denoted $\beta$ shows the gender-specific estimates of $\beta$ from
a probit in columns (1) and (2), and the gender-specific estimates of $\beta$
from our new method in columns (4) and (5). The estimated gap contributions of
each element of $\beta$, calculated via (11) on the full sample of individuals
surveyed, are given in column (3) for the uncorrected method and (6) for the
corrected method. We also report the effect of shifting all variables in $X$
for women to match their male means, as well as shifting all values of $\beta$
for women to match those of men. The “unexplained gap” is the additional
difference in estimated average aspirations among men and women that is not
explained purely by $X_{i}$ or $\beta$, which is partially driven by the
interaction of the gender gap in $X_{i}$ and the gender gap in $\beta$, and is
also partially driven by nonlinear effects of non-mean differences in the
distribution of $X_{i}$ between men and women.

In addition to the contributions of gender gaps in each variable and its
coefficient to the gap in entrepreneurial aspirations, we also report gender-
specific estimates of the correlation parameter, $\hat{\rho}$, in columns (4)
and (5), with the constrained values of zero implied by the uncorrected method
shown in columns (1) and (2). With multiple reminders, $\rho$ can be
estimated, so the p-value on $\hat{\rho}$ in the corrected model is a test of
the null hypothesis of no correlation in unobservables, as assumed by the
uncorrected method. This “overidentification-test” p-value is reported for men
and women in columns (1) and (2), with the overidentification-test p-value
from the test described in Section 2.4 shown in columns (4) and (5). We fail
to reject the null $\rho=0$ for both men and women at conventional
significance levels with p-values of 0.129 for both groups. We reject the
assumption of joint normality in (5) for men with a p-value of 0.040, while
failing to reject this assumption for women, with a p-value of 0.376.

Both methods estimate large gaps between men and women in entrepreneurial
aspirations, with the uncorrected and corrected gap estimates given by 0.207
and 0.194, respectively. We find marginally significant evidence of positive
selection bias for both men and women with p-values of 0.103 and 0.090,
respectively, which suggests that individuals who aspire to become
entrepreneurs may be more likely to respond to a survey about
entrepreneurship. The difference between the uncorrected and corrected gap
estimates of 0.013 (0.207-0.194) is also just shy of conventional significance
thresholds, with a p-value of 0.103, suggesting that the sample selection bias
might be larger for men than for women. This is consistent with a survey
response decision that is driven both by interest in entrepreneurship (higher
for men, lower for women) and conscientiousness (lower for men, higher for
women).

We do not find any particularly important predictors of entrepreneurial
aspirations among our predictor variables, with the total estimated
contribution of $X_{i}$ at 1.1 percentage points for both the uncorrected and
corrected models. The total uncorrected and corrected effects of $\beta$
(0.189 and 0.179) on the gender gap are overwhelmingly driven by the
contributions of the constants, which are 0.154 and 0.172 for the uncorrected
and corrected methods, respectively.

Overall, we find evidence of large gender gaps in entrepreneurial aspirations
among college students using a survey of undergraduate students at the
university of Wisconsin - Madison. The average gap among the survey
respondents is slightly smaller than the estimates provided by controlling for
observed determinants of survey response, with both of these slightly larger
than the gap estimates obtained using our method to control for otherwise
unobserved determinants of response using the timing of responses. We note
that we reject the bivariate normality assumption on model residuals for men,
so we take these results with a grain of salt. We plan to consider extensions
of our method in a future version of this paper that considers alternative
distributional assumptions when bivariate normality is rejected.

Table 2: Oaxaca-Blinder Decomposition of the Entrepreneurial Aspirations
Gender Gap

| Uncorrected | Corrected |
---|---|---|---
| Male | Female | Gap Portion | Male | Female | Gap Portion |
Variables | (1) | (2) | (3) | (4) | (5) | (6) |
Intention | 0.389 | 0.182 | 0.207 | 0.352 | 0.158 | 0.194 |
| (0.005) | (0.004) | (0.006) | (0.024) | (0.015) | (0.028) |
X |  |  |  |  |  |  |
Business Major | 0.128 | 0.096 | 0.004 | 0.121 | 0.085 | 0.004 |
| (0.002) | (0.002) | (0.000) | (0.001) | (0.001) | (0.000) |
STEM Major | 0.433 | 0.256 | -0.002 | 0.417 | 0.249 | -0.002 |
| (0.004) | (0.003) | (0.001) | (0.001) | (0.001) | (0.001) |
GPA | 3.408 | 3.438 | 0.005 | 3.288 | 3.392 | 0.004 |
| (0.004) | (0.003) | (0.001) | (0.002) | (0.001) | (0.001) |
ACT Math | 30.156 | 27.842 | -0.005 | 29.826 | 27.526 | -0.003 |
| (0.028) | (0.026) | (0.002) | (0.011) | (0.011) | (0.002) |
ACT Verbal | 54.864 | 55.810 | 0.002 | 54.026 | 54.819 | 0.002 |
| (0.086) | (0.072) | (0.000) | (0.033) | (0.030) | (0.000) |
Racial Minority | 0.268 | 0.243 | 0.001 | 0.301 | 0.271 | 0.001 |
| (0.003) | (0.003) | (0.000) | (0.001) | (0.001) | (0.000) |
International | 0.099 | 0.064 | 0.004 | 0.101 | 0.068 | 0.004 |
| (0.002) | (0.002) | (0.000) | (0.001) | (0.001) | (0.001) |
All X |  |  | 0.011 |  |  | 0.011 |
|  |  | (0.003) |  |  | (0.003) |
$\beta$ |  |  |  |  |  |  |
Business Major | 0.442 | 0.497 | -0.002 | 0.452 | 0.506 | -0.001 |
| (0.040) | (0.042) | (0.002) | (0.042) | (0.045) | (0.002) |
STEM Major | -0.033 | -0.057 | 0.001 | -0.040 | -0.047 | 0.000 |
| (0.028) | (0.032) | (0.003) | (0.031) | (0.035) | (0.003) |
GPA | -0.173 | -0.183 | 0.008 | -0.144 | -0.147 | 0.002 |
| (0.028) | (0.031) | (0.036) | (0.032) | (0.034) | (0.036) |
ACT Math | -0.005 | -0.009 | 0.028 | -0.005 | -0.006 | 0.000 |
| (0.004) | (0.004) | (0.044) | (0.004) | (0.005) | (0.040) |
ACT Verbal | -0.007 | -0.008 | 0.018 | -0.007 | -0.009 | 0.024 |
| (0.001) | (0.001) | (0.027) | (0.001) | (0.001) | (0.027) |
Racial Minority | 0.090 | 0.149 | -0.005 | 0.071 | 0.150 | -0.006 |
| (0.036) | (0.035) | (0.004) | (0.039) | (0.038) | (0.004) |
International | 0.290 | 0.543 | -0.006 | 0.344 | 0.505 | -0.004 |
| (0.055) | (0.060) | (0.002) | (0.061) | (0.065) | (0.002) |
Constant | 0.860 | 0.350 | 0.154 | 0.702 | 0.110 | 0.172 |
| (0.143) | (0.143) | (0.070) | (0.182) | (0.178) | (0.086) |
All $\beta$ |  |  | 0.189 |  |  | 0.179 |
|  |  | (0.007) |  |  | (0.028) |
All Unexplained Gaps |  |  | 0.006 |  |  | 0.004 |
|  |  | (0.005) |  |  | (0.005) |
Auxiliary Parameters |  |  |  |  |  |  |
$\rho$ | 0 | 0 | 0 | 0.061 | 0.064 | 0 |
| (.) | (.) | (.) | (0.040) | (0.042) | (.) |
Log-Likelihood | -299757.3 |  |  |  |  |  |
p-val overidentification | 0.129 | 0.129 |  | 0.040 | 0.376 |  |
p-val, Difference from Corrected | 0.103 | 0.090 | 0.103 |  |  |  |
Observations (responded) | 18227 | 22523 | 40750 | 18227 | 22523 | 40750 |
Observations (surveyed) | 124830 | 137524 | 262354 | 124830 | 137524 | 262354 |

_Notes:_ $X$ estimates in columns 1-2 are means among individuals who
responded to the survey, with columns 4-5 giving means among all individuals
surveyed. Columns 3 and 6 give uncorrected and selection-corrected estimated
contributions of each variable to the gap in entrepreneurial intentions from
equation (10). $\beta$ estimates give the estimated coefficients on variables
from a probit in columns 1-2, and from the selection correction in columns
4-5. Columns 3 and 6 give uncorrected and selection-corrected estimated
contributions of each value in $\beta$ to the gap in entrepreneurial
intentions from equation (11). p-values from the overidentification test
estimated via (7) are reported for men and women along with p-values on the
difference between corrected and uncorrected estimates of the average
entrepreneurial intention rate.

##  4 Conclusion

We develop a method for estimating population level relationships between
variables using survey data that corrects for selection into response that can
be used in settings with multiple reminders. The core insight is to estimate
differences in survey variables between early and late responders, and to
extrapolate from those differences to the rest of the population of interest,
for whom survey data does not exist. The method nests a standard Heckman
correction as a special case in which there is a single survey notification.
This means it corrects for both observed and unobserved determinants of
selection by permitting, but not relying on, the use of covariates that affect
both the variable surveyed and the survey response decision as well as those
that only affect the survey response decision, such as randomized response
incentives.

Our method is particularly noteworthy because of its low costs in multiple
dimensions. First, it does not require observation of any variables other than
those gathered during survey administration responses for any surveyed
individuals, setting it apart from methods that adjust for observable
determinants of survey response. Second, it does not require randomized survey
response incentives, which minimizes both administrative complexity and
financial costs, as well as making it ideal for secondary data from surveys
that were not administered with the goal of eliminating selection bias.
Finally, it relies on a straightforward alteration of an existing method that
is familiar to many applied researchers, implemented in statistical software
as a loop over Heckman corrections for each reminder.

Following Heckman (1979) and Van de Ven and Van Praag (1981), we model all
binary variables as arising from latent variable processes with normally
distributed residuals. The core of our method works with any assumed bivariate
distribution for residuals in the survey response decision and the outcome of
interest. The range of alternative methods includes methods that avoid
distributional assumptions at the cost of obtaining bounds rather than point
estimates on quantities of interest, such as that of Lee (2005) or Behaghel,
Crépon, Gurgand, and Le Barbanchon (2015). It also includes methods that allow
richer structure on the survey selection choice, such as the method of Dutz,
Huitfeldt, Lacouture, Mogstad, Torgovitsky, and Van Dijk (2021), at the cost
of requiring variables that affect the survey response decision but not the
outcome of interest in addition to survey reminders. We leave a survey of the
performance of alternative methods and distributional assumptions across a
variety of applications to future work, noting that the best method for a
particular application will be the one that leverages subject matter expertise
to obtain good in-sample fit and out-of-sample fit subject to constraints
imposed by data availability.

## References

  * (1)
  * Aldrich (2005) Aldrich, H. (2005): “Entrepreneurship,” _Research in the Sociology of Work_ , 15, 3–31.
  * Behaghel, Crépon, Gurgand, and Le Barbanchon (2015) Behaghel, L., B. Crépon, M. Gurgand, and T. Le Barbanchon (2015): “Please call again: Correcting nonresponse bias in treatment effect models,” _Review of Economics and Statistics_ , 97(5), 1070–1080.
  * Blinder (1973) Blinder, A. S. (1973): “Wage discrimination: reduced form and structural estimates,” _Journal of Human resources_ , pp. 436–455.
  * Canning, Haque, and Wang (2012) Canning, J., M. Haque, and Y. Wang (2012): “Women at the wheel: Do female executives drive start-up success,” _Dow Jones and Company_.
  * DiNardo, Matsudaira, McCrary, and Sanbonmatsu (2021) DiNardo, J., J. Matsudaira, J. McCrary, and L. Sanbonmatsu (2021): “A Practical Proactive Proposal for Dealing with Attrition: Alternative Approaches and an Empirical Example,” _Journal of Labor Economics_ , 39(S2), S507–S541.
  * Dutz, Huitfeldt, Lacouture, Mogstad, Torgovitsky, and Van Dijk (2021) Dutz, D., I. Huitfeldt, S. Lacouture, M. Mogstad, A. Torgovitsky, and W. Van Dijk (2021): “Selection in Surveys: Using Randomized Incentives to Detect and Account for Nonresponse Bias,” Discussion paper, National Bureau of Economic Research.
  * Fairlie (2005) Fairlie, R. W. (2005): “An extension of the Blinder-Oaxaca decomposition technique to logit and probit models,” _Journal of economic and social measurement_ , 30(4), 305–316.
  * Fisher (1922) Fisher, R. A. (1922): “On the mathematical foundations of theoretical statistics,” _Philosophical transactions of the Royal Society of London. Series A, containing papers of a mathematical or physical character_ , 222(594-604), 309–368.
  * Gneezy and Rustichini (2000) Gneezy, U., and A. Rustichini (2000): “Pay enough or don’t pay at all,” _The Quarterly journal of economics_ , 115(3), 791–810.
  * Greene, Hart, Gatewood, Brush, and Carter (2003) Greene, P. G., M. M. Hart, E. J. Gatewood, C. G. Brush, and N. M. Carter (2003): “Women entrepreneurs: Moving front and center: An overview of research and theory,” _Coleman White Paper Series_ , 3(1), 1–47.
  * Guzman and Kacperczyk (2019) Guzman, J., and A. O. Kacperczyk (2019): “Gender gap in entrepreneurship,” _Research Policy_ , 48(7), 1666–1680.
  * Heckman (1979) Heckman, J. J. (1979): “Sample selection bias as a specification error,” _Econometrica: Journal of the econometric society_ , pp. 153–161.
  * Lee (2005) Lee, D. S. (2005): “Training, wages, and sample selection: Estimating sharp bounds on treatment effects,” .
  * Molina Millán and Macours (2017) Molina Millán, T., and K. Macours (2017): “Attrition in randomized control trials: Using tracking information to correct bias,” .
  * Oaxaca (1973) Oaxaca, R. (1973): “Male-female wage differentials in urban labor markets,” _International economic review_ , pp. 693–709.
  * Van de Ven and Van Praag (1981) Van de Ven, W. P., and B. M. Van Praag (1981): “The demand for deductibles in private health insurance: A probit model with sample selection,” _Journal of econometrics_ , 17(2), 229–252.