One method to make waveguide or fiber couplers is to use straight sections of the guides where the evanescent modes of one guide overlap with the modes of a second guide, eg, a directional coupler. The light from one guide slowly transfers back and forth between the guides. By choosing the distance between the guides and the length of the coupling region, it is possible to couple any desired fraction of the light from one waveguide to the other. The coupling strength can be very sensitive to the distance between the guides and it is important to ultimately choose a design that can function acceptably given the types of imperfections that are expected for your manufacturing process.

The file waveguide_coupler.lms has two Silicon on Insulator (SOI) waveguides that are 500nm wide and 200nm thick. The space between the waveguides is 50nm. There are two mesh override regions used:
1.One region covers both guides and extends about 100nm outside them, that reduces the size of the mesh by a factor of 4 in order to accurately resolve the modes inside the guides and in the evanescent field; and
2.A second region that sets the values of dy between the guides to 2nm. This is important because the coupling of the guides can be quite sensitive to the space between them.
We can first calculate the coupling length by considering the difference in index between the 2 coupled modes, as shown below.
The index difference is
$$ \Delta n=n_{1}n_{2}=2.3780222.317864=0.0601584 $$
and, if we initially have 100% of the power in waveguide 1, then the power in waveguide 2 is given by
$$ P_{2}(L)=P_{0} \sin ^{2}\left(\frac{\pi L \Delta n}{\lambda_{0}}\right) $$
where L is the length of the coupling region, P0 is the input power, and λ0 is the free space wavelength. It is easy, for example, to calculate the coupling length to achieve any desired power coupling to waveguide 2 by solving for L:
$$ L=\frac{\lambda_{0}}{\pi \Delta n} \sin ^{1}(\sqrt{\frac{P_{2}}{P_{0}}}) $$
For example, the length required to couple 100% of the light from waveguide 1 to waveguide 2 is simply
$$ L_{100 \%}=\frac{\lambda_{0}}{2 \Delta n}=\frac{1.55 \mu m}{2 \cdot 0.0601584}=12.8827 \mu m $$
How does this work?
Consider the situation where light is propagating in waveguide 1 and we introduce a second waveguide abruptly, as shown in the left image below. The Ey component of the mode propagating in waveguide 1 initially is shown in the plot on the right. The effective index of the initial mode is about n0 = 2.323.
Once the second waveguide is introduced, we now have 2 modes, which are approximately given by the sum and difference of the modes of the original waveguides, and have effective indices n1 = n0 + Δn/2 and n2 = n0  Δn/2. If we consider the sum of these modes we have
so we see that the original waveguide mode matches well with the sum of the 2 modes that exist when both waveguides are present. The original waveguide mode therefore excites both of these modes equally, but they propagate with different phase velocities. The E field at a distance L from the point where the second waveguide is introduced is given by
$$ \begin{aligned} \vec{E}(L) &=\vec{E}_{1} \exp \left(\frac{2 \pi i n_{1}}{\lambda_{0}} L\right)+\vec{E}_{2} \exp \left(\frac{2 \pi i n_{2}}{\lambda_{0}} L\right) \\ & \approx\left(\vec{E}_{1}^{0}+\vec{E}_{2}^{0}\right) \exp \left(\frac{2 \pi i n_{1}}{\lambda_{0}} L\right)+\left(\vec{E}_{1}^{0}\vec{E}_{2}^{0}\right) \exp \left(\frac{2 \pi i n_{2}}{\lambda_{0}} L\right) \\ &=2 \vec{E}_{1}^{0} \exp \left(\frac{2 \pi i n_{0}}{\lambda_{0}} L\right) \cos \left(\frac{\pi \Delta n}{\lambda_{0}} L\right)+2 i \vec{E}_{2}^{0} \exp \left(\frac{2 \pi i n_{0}}{\lambda_{0}} L\right) \sin \left(\frac{\pi \Delta n}{\lambda_{0}} L\right) \end{aligned} $$
where E1 and E2 are the modes of the coupled waveguide system, and E10 and E20 are the unperturbed modes of the single waveguide positioned at location 1 or location 2 respectively. We see that the light propagates like the unperturbed waveguide system, but the oscillates back and forth between the 2 guides periodically due to the sine and cosine terms.
In reality, we rarely introduce a second waveguide abruptly and we generally don't remove it abruptly. A real device might look something like this
In many cases, the effect of the bent region can be ignored. The effect of the bent region can be investigated in detail either using the propagator in MODE, or if the waveguide is small enough, it may be possible to run a 3D FDTD simulation.
Reconstructing an image of the propagating field
It may be useful to reconstruct the image of the field as it propagates over large distances. The script file waveguide_couplers.lsf can be used to do this after opening the file waveguide_couplers.lms. It first deletes the right waveguide and calculates the field of the input guide which is copied to a global dcard called "E0", and looks like the mode below.
The right waveguide is then reintroduced by the script and the modes of the coupled waveguide system are calculated. The mode "E0" is then decomposed onto the 2 modes of the coupled system and propagated an arbitrary distance using the propagate command. This is repeated for 100 different lengths from 0 to 40 microns, and the final figure is created which shows the electric field intensity vs x and L at y=0.1 microns.
We can see that the length for 100% coupling is about 13 microns, as expected.
It is also possible to simulate the same structure using the propagator. The simulation file waveguide_coupler.lms also contains a propagator region which can be activated by right clicking on the propagator region on the object tree.
The propagator simulation contains a mode source. If you choose to user select a mode, you will get the mode solver shown below. As shown in the eigenmode solver example above, the propagation length is related to the effective index.
From the screenshot below you can see that the index difference is
$$ \Delta n=n_{1}n_{2}=2.4754232.391182=0.084241 $$
Notice that this effective index difference is slightly higher than when we used the eigenmode solver, therefore running the varFDTD as is will return a coupling length noticeably shorter than expected. Clearly the FDE solver gives a more accurate effective index difference for this structure because it calculates full 2D mode profiles. In contrast, the propagator first compresses the z dimension and then calculates the modes using the compressed data.
It is very important to return the correct effective index as it directly affects the coupling length. In varFDTD simulation, we can slightly modify the spacing between the two waveguides in order to get the correct effective index difference. A index difference of ~ 0.06 can be attained by setting the waveguide spacing to 70nm, using the code below.
# set the 70 nm gap
gap = 70e9;
setnamed("left", "y", ( 500e9 + gap )/2);
setnamed("right", "y", ( 500e9 + gap )/2);
setnamed("mesh_gap", "y span", gap);
Then, we can set the x max of the right waveguide to 20 microns, and set the integrated mode source so that it selects the fundamental mode. If we run the simulation, and send the monitor1 data to the frequency monitor, we can see that the coupling length is approximately 12 microns.
# introducing the abrupt waveguide and run the simulation
setnamed("right", "x max", 20e6);
run;
The oscillations seen in the electric field intensity above are due to the fact that there is an additional mode excited. The Effective Index tab for the propagator contains an option to clamp the effective index values to the physical material properties. The extra mode which is excited in the image above can be removed by unchecking this option.
Since removing the material clamping changes the variational effective index method, the difference between the real part of the effective indices of the modes changes again. Therefore, in order to obtain the image below, the mode solver was used to find the correct spacing between the waveguides. After removing the clamping a spacing of 100nm gave a refractive index difference of 0.06. If we run the simulation, and send the monitor1 data to the frequency monitor, we can see that the coupling length is approximately 12.5 microns.
# disabling the clamp index option and set the 100 nm gap
setnamed("varFDTD", "clamp values to physical material properties", 0);
gap = 100e9;
setnamed("left", "y", ( 500e9 + gap )/2);
setnamed("right", "y", ( 500e9 + gap )/2);
setnamed("mesh_gap", "y span", gap);
# introducing the abrupt waveguide and run the simulation
setnamed("right", "x max", 20e6);
run;
It is also possible to simulate the same structure using the EME. The simulation file waveguide_coupler.lms also contains a EME region which can be activated by right clicking on the appropriate solver region on the object tree.
In EME, we can use port object as a mode source, find the two coupled modes and calculate the coupling length. Nevertheless, for the simulation itself, we need to inject fundamental TE mode into one waveguide only. To do this, setup the "x min" position of the right waveguide to 24.5um. This will allow us to use the port offset settings as shown on the figure below to calculate the injected mode in the region between x=25 and 24.5um where we have only the left waveguide now.
The EME results extracted from the field monitor show that the simulated coupling length is 12.7um, which is in good agreement with the other solvers and analytical calculation.
We have demonstrated above that we can use three different solvers for the evanescent coupler simulation and get results that agree with analytical calculation. Hence, we might be interested in simulation speed comparison to choose the solver that provides us the results fastest for this specific simulation using a modest desktop machine.
Solver 
Simulation Time 
Coupling Length 
FDE and analytical calculations 
N/A 
12.88 um 
FDE and propagation 
~20 s 
13 um 
EME (run+propagate) 
~10 s 
12.7 um 
varFDTD 
~60 s 
12.5 um * 
* This varFDTD coupling length results is obtained after the adjustment on the gap distance to achieve the same delta_neff, as explained in the above section.
Note: Accuracy The eigenmode solvers should in principle give a more accurate result than varFDTD because they calculate the full 2D mode profiles, whereas the varFDTD solver is less accurate since it neglects the coupling in the z dimension. 
For small variations, it should be still possible to use the FDE to calculate the neff and modes, and get a rough idea of the coupling length. However, the propagation command assumes uniform waveguide crosssection, and therefore it does not account for any nonuniform regions. To fully account for the variation in the coupling region, it will be necessary to use some propagation methods simulate the whole device, such as EME, varFDTD, or FDTD.