The Hong–Ou–Mandel effect can be described like this: Two photons enter a particular type of beamsplitter, one from each side. If the two photons are "identical" in the sense defined below, then the probability of detecting one photon on each side of the output is zero. Both photons will always be detected on the same side of the output, though we cannot predict which side.
The notation used in the OP is partly unclear to me, but one thing is clear: In a state with two photons, shifting one photon's phase by a constant $\theta$ cannot affect anything. It merely multiplies the two-photon state-vector by and overall factor of $\exp(i\theta)$. Physical predictions depend only on the ray (the one-dimensional Hilbert space spanned by the given state-vector), not on the individual state-vector, so multiplying the state-vector by $e^{i\theta}$ has no physical effect. (This is in contrast to the effect of a relative phase between two terms in a superposition.) This is restated below in mathematical notation.
The OP considered only plane waves, but the following calculation allows the two incoming photons to have arbitrary longitudinal profiles, described by complex-valued functions $f$ and $g$, respectively. The HOM effect occurs whenever $f\propto g$, which includes the case $f=e^{i\theta}g$ considered in the OP. I'll use the creation/annihilation operator formalism, because this makes the calculation easy and make the reason for the result very clear.
Calculation
Let $a^\dagger(k)$ denote the operator that, when applied to any state-vector, adds a photon traveling diagonally downward with wavenumber $k$.
Let $b^\dagger(k)$ denote the operator that, when applied to any state-vector, adds a photon traveling diagonally upward with wavenumber $k$.
Let $|0\rangle$ denote the state with no photons (vacuum state). This state is annihilated by $a(k)$ and $b(k)$, the adjoints of $a^\dagger(k)$ and $b^\dagger(k)$.
The definitions of $a^\dagger$ and $b^\dagger$ are illustrated here:
The letter ($a$ or $b$) indicates the photon's direction, and the argument $k$ indicates its wavenumber along that direction. For any complex-valued function $f$, use the abbreviation
$$
a^\dagger(f) := \int dx\ f(k)a^\dagger(k)
\tag{1}
$$
and similarly for $b^\dagger(f)$. The commutation relations are
\begin{gather}
[a(f),a(g)]
= [b(f),b(g)] = 0
\\
[a(f),b(g)]
=[a(f),b^\dagger(g)] = 0
\tag{2}
\end{gather}
and
$$
[a(f),a^\dagger(g)] =
[b(f),b^\dagger(g)] = \int dk\ f^*(k) g(k).
\tag{3}
$$
The effect of an ideal beamsplitter can be modeled using the transformation
\begin{align}
a^\dagger(f) &\to \big(a^\dagger(f)+ib^\dagger(f)\big)/\sqrt{2}
\\
b^\dagger(f) &\to \big(b^\dagger(f)+ia^\dagger(f)\big)/\sqrt{2}
\tag{4}
\\
|0\rangle &\to |0\rangle.
\end{align}
This transformation has these properties:
It is unitary.
It is symmetric in $a$ and $b$. In other words, the beamsplitter has up/down symmetry.
It splits the photon's direction but otherwise preserves the wavenumber $k$.
The factors of $i$ are necessary in order for the transformation to be both unitary and symmetric in $a$ and $b$. That's important, because the HOM effect relies on these factors of $i$.
(The WP article linked in the OP uses a different unitary transformation that doesn't have factors of $i$ and is not manifestly symmetric in $a$ and $b$. That transformation can be converted into this one by absorbing a factor of $i$ into either $a$ or $b$, which doesn't affect the commutation relations. And as ZeroTheHero indicated in a comment, (4) is not the most general form for a 50:50 beamsplitter, but it is sufficient for the present purposes.)
To derive the HOM effect, suppose that the state prior to the beamsplitter is
$$
a^\dagger(f) b^\dagger(g)|0\rangle.
\tag{5}
$$
This is a two-photon state, with one photon traveling diagonaly downward and one traveling diagonally upward. Since this is the initial state, the photons are understood to be approaching the beamsplitter. (We could make this explicit by using the machinery of quantum field theory — local field operators — to relate the wavenumber domain to the position domain.) According to the transformation defined above, the effect of the beamsplitter is
\begin{align}
a^\dagger(f) b^\dagger(g)|0\rangle
&\to
\frac{1}{2}
\big(a^\dagger(f)+ib^\dagger(f)\big)
\big(b^\dagger(g)+ia^\dagger(g)\big)|0\rangle
\\
&=
\frac{1}{2}
\big(
ia^\dagger(f)a^\dagger(g)
+ib^\dagger(f)b^\dagger(g)
\\
&\phantom{\frac{1}{2}\big(}
+a^\dagger(f)b^\dagger(g)
- a^\dagger(g)b^\dagger(f)
\big)|0\rangle.
\tag{6}
\end{align}
This final state has four terms:
The $aa$ and $bb$ terms, in which both photons exit the beamsplitter in the same direction,
The $ab$ and $ba$ terms, in which the two photons exit the beamsplitter in different directions.
Conclusion
In the special case $f\propto g$, a glance at the result (6) shows that the $ab$ and $ba$ terms cancel each other. This is the HOM effect. The relative minus sign between these two terms comes from $i^2=-1$. The linked WP article words it this way:
The Hong–Ou–Mandel effect ... occurs when two identical single-photon waves enter a 1:1 beam splitter, one in each input port. When the photons are identical, they will extinguish each other. If they become more distinguishable, the probability of detection will increase.
That wording is unclear. Here's the decoder ring:
"Identical" means $f\propto g$. "Distinguishable" means that $f$ and $g$ are not proportional to each other. This wording apparently alludes to the single-photon states $a^\dagger(f)|0\rangle$ and $a^\dagger(g)|0\rangle$, which are physically distinguishable if and only if $f$ and $g$ are not proportional to each other.
The "probability of detection" apparently means the probability of detecting photons at both output ports, which cannot happen when $f\propto g$. (That's what the preceding calculation showed.) This wording is probably related to this statement from https://arxiv.org/abs/1711.00080: "In HOM interference, we are often interested in the coincidence probability, that is, the probability of detecting one photon in each output port of the beam splitter."
The OP considers the case $f=e^{i\theta}g\propto g$. The preceding calculation shows that the $ab$ and $ba$ terms cancel each other whenever $f\propto g$. This agrees with the conclusion of the OP's calculation.
What does "coherent" mean?
The most important point is that the result cannot depend on $\theta$, and this is already evident from the initial state (5), without calculating anything at all. The replacement $g\to e^{i\theta}g$ can't have any effect, because it merely multiplies the initial state-vector by an overall factor:
$$
a^\dagger(f) b^\dagger(g)|0\rangle
\to
e^{i\theta}a^\dagger(f) b^\dagger(g)|0\rangle.
\tag{7}
$$
Physical predictions depend only on the ray (the one-dimensional Hilbert space spanned by the given state-vector), not on the individual state-vector.
Contrast this to the single-photon state with a $\theta$-dependent relative phase:
$$
a^\dagger(f) |0\rangle +
e^{i\theta} b^\dagger(g)|0\rangle.
\tag{8}
$$
In this case, the value of $\theta$ does matter. In the OP, the phase shift $g\to e^{i\theta}g$ is described as making the two photons "not [mutually] coherent." Not sure why such a word would be used in the context of (7), but the word does make sense in the context of a state like
$$
\exp\big(b^\dagger(g)\big)|0\rangle
=\sum_{n\geq 0} \frac{\big(b^\dagger(g)\big)^n}{n!}|0\rangle,
\tag{9}
$$
which is often used as a model of the light emitted by a laser. This so-called "coherent state" is affected by the phase shift $g\to e^{i\theta} g$. If we consider a situation in which the input to the beamsplitter consists of two of these "laser beams," then the wording used in the OP would make more sense. Explicitly, the input state in that case would be
$$
\exp\big(a^\dagger(f)\big)
\exp\big(e^{i\theta}b^\dagger(g)\big)|0\rangle,
\tag{10}
$$
which is equivalent to
$$
\exp\big(a^\dagger(f)+e^{i\theta}b^\dagger(g)\big)|0\rangle.
\tag{11}
$$
In contrast to the two-photon state (7), the phase $\theta$ does affect the physical significance of the state (10)-(11), as it does in (8).
